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Resolution  of  Source  Functions  and  Structure  from  Near-field  Data 


J.  S.  Barker,  R.  W.  Burger  and  L.  J.  Burdick 
Woodward-Clyde  Consultants,  566  El  Dorado  St.,  Pasadena,  CA  91101 

INTRODUCTION 

An  important  problem  in  nuclear  explosion  seismology  is  the  accurate 
determination  of  the  reduced  displacement  potential  (RDP).  This  paper 
summarizes  the  results  of  three  studies  concerning  the  determination  of  the 
RDP  from  near-field  data,  and  implications  for  teleseismic  attenuation 
estimates. 

THE  FORM  OF  THE  RDP 

Two  representations  of  the  RDP  are  currently  in  use,  differing  principally 
in  their  high-frequency  asymptotic  behavior  and  in  their  scaling  with  yield. 
The  displacement  spectra  derived  from  these  RDPs  for  the  Pahute  Mesa  event 
HALFBEAK  are  shown  in  Figure  1.  At  high  frequencies,  the  Mueller-Murphy 
(1971,  M-M)  spectrum  falls  off  as  f  ,  whereas  the  Helmbergcr-Hadley  (1981, 

__  O 

H-H)  spectrum  falls  off  as  f  .  In  the  1-5  Hz  band  resolvable  with  near-field 
data,  however,  the  two  spectra  are  virtually  indistinguishable.  This  is 
perhaps  to  be  expected,  since  both  formalisms  were  derived  by  fitting 
near-field  data  within  this  band.  A  comparison  of  the  spectral  slopes  and  the 
spectral  ratio  of  these  source  functions  (Figure  1)  indicates  that  they 
diverge  significantly  only  above  about  5  Hz  and  below  1  Hz. 

Comparisons  of  elastic  generalized  ray  synthetic  seismograms  computed  for 
both  RDP  representations  with  observed  near-field  data  are  given  in  Barker  et 
al  (1986),  but  a  typical  example  is  shown  in  Figure  2.  Shown  at  the  top  is 
the  observed  vertical  velocity  waveform  from  the  event  BOXCAR  recorded  at  a 
range  of  10  km.  Shown  in  the  dashed  line  is  the  H-H  synthetic  using 
parameters  determined  by  Barker  et  al  (1985),  and  in  the  dotted  line  is  the 
M-M  synthetic  scaled  to  the  announced  yield  of  1300  kt.  The  first  peak 
amplitude  and  frequency  content  of  the  M-M  synthetic  are  too  high  to  match  the 
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observations.  Shown  below  the  time-domain  velocity  waveforms  are  the 
amplitude  spectra  for  the  window  indicated  above  the  traces.  Only  those 
portions  of  the  spectra  from  1-10  Hz  have  been  displayed.  Once  again,  from 
1-5  Hz,  the  two  synthetic  spectra  are  virtually  identical  and  differ  from  the 
observed  spectrum  mostly  due  to  later,  unmodeled  arrivals  in  the  window. 

Above  5  Hz,  the  M-M  spectrum  is  enhanced  relative  to  the  observations,  while 
the  H-H  spectrum  continues  to  track  the  data. 

Since  the  elastic  M-M  synthetic  is  too  rich  in  high  frequency  energy  to 
satisfy  the  observations,  we  investigated  whether  the  effects  of  local 
anelastic  attenuation  could  improve  the  fit.  Figure  3  shows  a  comparison  in 
the  time  and  frequency  domains  between  elastic  H-H  and  M-M  synthetics  and  M-M 
convolved  with  a  Futterman  (1962)  average  attenuation  operator  for  a  Q  of  100 
(MMQ) .  The  synthetics  are  computed  for  a  1000  kt  source  at  10  km,  similar  to 
the  BOXCAR  comparison  shown  in  Figure  2.  The  value  of  Q=100  brings  the  first 
peak  amplitude  of  the  MMQ  synthetic  into  agreement  with  the  H-H  amplitude,  and 
has  the  desired  effect  of  depressing  the  high-frequency  spectral  amplitudes. 
Barker  et  al  (1986)  indicate  that  Q  values  as  high  as  100  are  too  high  for  the 
M-M  synthetics  to  match  observations  at  shorter  ranges.  We  conclude  that  with 
a  moderate  level  of  local  anelastic  attenuation,  the  M-M  source  representation 
is  indistinguishable  from  the  elastic  H-H  source,  but  that  a  range- dependent, 
or  more  precisely,  a  depth-dependent  attenuation  model  must  be  specified 
before  the  M-M  synthetics  can  be  used  with  near-field  data. 

IMPLICATIONS  FOR  TELESEISMIC  ATTENUATION  ESTIMATES 

It  is  apparent  from  Figure  1  that  the  spectral  amplitudes  of  the  H-H  and 
M-M  sources  are  significantly  different  in  the  0.5-2  Hz  band  in  which 
teleseismic,  short-period  P  waves  are  observed.  Within  this  band,  the  scaling 
relations  result  in  the  M-M  source  having  higher  absolute  amplitudes  than  the 
H-H  source.  Thus,  the  former  requires  higher  t  values  to  match  observed 
teleseismic  amplitudes.  The  figure  also  shows  that  out  to  about  4  Hz,  the  two 
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spectra  have  comparable  spectral  roll-off.  The  more  rapid  decay  of  the  H-H 
source  occurs  only  above  4  Hz.  We  have  used  the  near-field  results  to  model 
the  amplitudes  and  periods  of  teleseismic,  short-period  waveforms  in  order  to 
determine  estimates  of  teleseismic  attenuation.  Both  frequency- independent 
(Futterman,  1962)  and  frequency-dependent  (Minster,  1978)  models  were 
considered.  Details  of  the  study  may  be  found  in  Burger  et  al  (1986),  but  an 
overview  is  given  in  Figure  4.  Shown  are  various  source  amplitude  spectra 

9 

(multiplied  by  f  )  convolved  with  attenuation  operators.  The  top  curve  is  the 

model  inferred  by  Der  et  al  (1985)  based  on  matching  the  teleseismic  spectral 

* 

shape  in  the  0.5-4  Hz  band  using  a  VSB  source  and  at  of  0.34  s.  The  other 

curves  result  from  short-period  P-wave  modeling,  and  include  the  H-H  source 

* 

with  a  Minster  Q  operator  with  t  =1.0  s  and  TAUm=0.044  s  (H  M) ,  the  M-M  source 
with  Minster  Q  parameters  t*=1.25  s  and  TAUm=0.051  s  (MM),  the  H-H  source 
with  a  Futterman  t*  of  0.79  s  (H  F)  and  the  M-M  source  with  Futterman  t  of 
0.99  s  (M  F).  It  is  immediately  apparent  that  the  Der  et  al  (1985)  model  does 
not  fit  the  time-domain  amplitudes  near  1  Hz.  Similarly,  both  the  H-H  and  M-M 
sources  with  the  Futterman  attenuation  operator  have  spectral  fall-off 
inconsistent  with  the  data  modeled  by  Der  et  al  (1985).  However,  both  of 
these  source  models  with  the  frequency-dependent  Minster  Q  operator  match  both 
the  time-domain  amplitudes  and  spectral  fall-off,  within  a  reasonable  scatter 
in  the  data. 

THE  EFFECT  OF  LOCAL  STRUCTURE 

The  determination  of  source  parameters  from  near-field  data  is  sensitive  to 
errors  or  inadequacies  in  the  assumed  local  velocity  structure  model.  For 
example.  Stump  and  Johnson  (1984)  obtained  time-varying  moment  tensors  for 
three  events,  interpreting  the  deviatoric  parts  in  terms  of  non-isotropic 
source  radiation.  These  may  just  as  well  have  resulted  from  purely  isotropic 
sources  but  inaccurate  Green's  functions.  In  order  to  quantify  and  account 
for  these  tradeoffs,  we  have  developed  a  simultaneous  inversion  for  the 
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parameters  of  the  explosion  source  and  the  velocity  structure  models.  Figure 
5  shows  the  results  of  a  test  inversion  using  a  synthetic  data  set.  The 
synthetic  data  were  generated  for  a  2-gradient  model  approximating  the  Pahute 
Mesa  model  of  Hartzell  et  al  (1983).  In  the  inversion,  the  structure 
parameters  are  the  gradient,  the  velocity  at  the  top  of  the  layer,  and  the 
depth  to  the  top  of  the  layer.  The  gradient  is  discretized  into  a 
plane-layered  model  based  on  a  critical  wavelength,  and  synthetics  are 
computed  using  asymptotic  generalized  ray  theory.  For  this  case,  the  source 
is  represented  by  the  parameters  of  the  H-H  source  (holding  the  B  parameter 
constant  at  1.0).  The  residuals  are  the  normalized  cross-correlation 
coefficient  (a  measure  of  waveform  fit),  the  difference  of  normalization 
factors  (a  measure  of  amplitude)  and  the  time  lag  to  the  peak  of  correlation 
(a  measure  of  travel  time).  Partial  derivatives  are  computed  numerically  by 
perturbing  each  of  the  starting  parameters,  and  a  linearized  generalized 
inverse  is  solved.  The  results  are  listed  in  Table  1 .  For  the  first  five 
iterations,  only  the  parameters  of  the  source  and  upper  gradient  layer  were 
allowed  to  vary.  The  waveforms  in  Figure  5  indicate  rapid  convergence  to  the 
upper  gradient,  but  details  of  waveforms  at  farther  stations  are  not  well 
modeled,  and  the  source  parameters  are  still  unresolved.  The  parameters  of 
the  lower  gradient  could  not  be  resolved  using  only  stations  out  to  12  km. 
However,  when  artificial  data  out  to  20  km  was  included,  these  parameters 
approach  their  correct  value 3  in  iterations  6-10.  The  parameter  resolution 
matrix  from  the  inversion  indicates  the  extent  to  which  source  and  structure 
parameters  trade-off  as  a  function  of  iteration.  In  addition,  with  an 
estimate  of  the  variance  of  the  residuals  using  real  data,  error  ellipsoids 
for  the  inversion  parameters  may  be  formed  and  used  to  evaluate  the  quality  of 
the  solution. 
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Figure  1.  (Left)  Predicted  source  displacement  spectra  for 
HALFBEAK.  The  Mue 1 1 e r - Mur p hy  source  models  are 

based  on  the  announced  yield.  The  Helmberger-Hadley  are  based 
on  the  near  field  modeling  results  of  Barker  et  al.  (1985). 
(Center)  Spectral  slopes  of  the  displacement  spectra  versus 
frequency.  (Right)  The  spectral  ratio  of  the  two  source 
mode  1 s . 
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Figure  4.  A  comparison  of  the  theoretical  spectra 
for  the  Helmberger-Hadley  and  Mueller-Murphy  sources 
convolved  with  the  Futterman  and  Minster  Q  operators 
to  the  spectrum  preferred  by  Der  et  al.  (1986)  for 
Pahute  Mesa  (von  Seggern  and  Blandford  source  with 
Futterman  Q  operator,  t*=0.34). 
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Figure  5.  Inversion  waveforms  and  structure  model  results.  The 
starting  model  and  results  of  selected  iterations  are  shown  and 
compared  with  the  synthetic  data  set.  For  the  first  5  iterations 
only  receivers  to  12  km  were  included.  For  iterations  6-10,  the 
inclusion  of  receivers  to  20  km  allowed  the  lower  gradient  to  be 
resolved.  Arrows  indicate  the  time  windows,  and  peak  amplitudes 
are  shown . 
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Introduction 

The  method  used  in  that  study  is  derived  from  the  discret  waves 
numbers  method.lt  has  been  used  to  characterize  the  effects  of 
different  large  scale  heterogeneities  on  the  Lg  waves  propagation. 

Two  types  of  effect  have  been  studied. 

First,  the  effects  due  to  Moho  irregularities.  The  Lg  waves  train 
is  a  superposition  of  reflected  S  waves  trapped  within  the  crust. 
Consequently  any  variation  of  the  Moho  depth  might  affect  the  waves 
guide  quality  and  leads  to  some  apparent  strong  attenuation  zone  for 
Lg  waves . 

Second,  we  look  at  the  effects  of  sedimentary  basins  on  the  Lg  waves 
propagation  path. 

The  study  is  limited  to  the  SH  in  two  dimensions.  Irregularities 
dimensions  are  selected  to  fit  the  situations  met  in  central  part  of 
France  for  which  general  attenuation  features  have  been  studied 
(  Campillo  ,  Plantet  and  Bouchon  1985). 


10 


I  Variation  of  the  Moho  depth 


The  crust  model  is  given  in  Figure  1.  Computations  have  been  made 
for  seismic  sources  at  10  km  depth,  located  at  0,  50,  100,  150  and  200  km 
from  the  Moho  irregularity. 

Figure  2  shows  maximum  amplitude  variations  for  Lg  waves  versus 
distance,  for  a  planar  crust  mode  (  thick  solid  line)  and  for  a  crust 
model  with  irregularity  at  these  previous  distances. 

In  the  irregularity  vicinities,  steepen  reflectors  bring  strong  amplitude 
variations . 

On  the  contrary,  at  larger  distances,  no  important  and  systematic 
differences  are  seen  between  a  Moho  without  or  with  irregularity. 

This  points  out  that  any  Moho  ascent  should  commonly  bring  limited 
attenuation  effects  on  Lg  waves. 

II  Sedimentary  basin 

We  look  now  at  the  influence  of  a  thicken  sedimentary  layer  or  a  sedi¬ 
mentary  basin  as  described  in  Figure  3  on  Lg  waves  at  the  same  distances 
from  the  irregularity(as  describe  in  Figure  1  )  seismic  sources  are  always 
at  10  km  depth. 

The  synthetic  seismograms  are  presented  in  Figure  4  for  various 
distances  from  0  to  290  km,  every  10  km,  for  both  a  planar  crust  (  upper 
seismograms  )  and  a  thicken  sedimentary  layer  at  the  top  of  the  crust 
(  lower  seismograms). 

Clearly  seen  is  the  formation  of  a  Love  wave  in  the  case  of  a  sedimentary 
basin.  As  superficial  layers  are  propagation  structures  which  strongly 
affect  Love  waves.,  the  converted  Lg  waves  to  Love  waves  energy  will  be 
rapidly  dissipated. 

Within  the  sedimentary  basin  itself,  classical  amplification  phenomena 
might  happen. 

These  two  rough  synthetic  examples  suggest  to  verify  their  conclusions 
for  some  real  cases. 
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Conclusion 


These  results  pointed  out  the  importance  of  graben  zones  for  Lg  waves 
propagation.  Nevertheless,  at  this  stage  of  the  study  we  do  not  possess 
any  separation  tool  between  two  types  of  phenomena  which  could  create 
attenuation  : 

-  small  scale  diffraction  (  or  scattering  ) 

-  large  scale  geometrical  effect. 

A  study  of  frequency  dependancy  for  Q  should  help  in  the  interpretation 
of  our  data. 

On  the  other  hand,  developments  of  numerical  simulation  methods 
encourage  us  to  investigate  more  complex  zones  in  order  to  interprete 
the  quasi-extinction  phenomena  of  Lg  waves  we  have  observed  along  peculiar 
paths. 
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Figure  1 


crust  model  with  Moho  irregularity  -  the  star  is  the 


Figure  2 


Figure  3 


Figure  4 


Figure  5 


Figure  6 


Figure  7 


seismic  source  -  seismograms  are  built  at  various  distances, 
maximum  amplitudes  of  Lg  waves  versus  distance  (  2D  ). 
crust  model  with  a  sedimentary  basin 

synthetics  obtained  with  a  planar  crust  model  (  high)  and  with 
a  sedimentary  irregularity  (  low  ) 

Rays  paths  used  to  build  up  a  map  of  Q  G  factor  for  France 
Map  of  confidence  index 

Map  of  France  (  Center  )  for  Q  factor  at  1  Hz 
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AFGL/DARPA  REVIEW  OF  NUCLEAR  TEST  MONITORING  BASIC  RESEARCH 
U.S.  AIR  FORCE  ACADEMY,  6-8  MAY  1986 

Modeling  Explosions  Using  2-D  Numerical  Methods 

Donald  V.  Helmberger,  Richard  Stead  and  John  Vidale 
Seismological  Laboratory 
California  Institute  of  Technology 

Introduction 

It  is  well  known  that  yield  versus  m^  plots  often  show 
geographical  correlations,  for  example  see  Lay  et.  al  (1984). 
Some  of  this  variation  could  be  caused  by  depth  effects  and, 
perhaps,  tectonic  release.  Another  possibility  is  scattered 
Rayleigh  waves.  Alewine  has  shown  that  many  shots  near  the  edge 
of  Yucca  Flats  and  Pahute  Caldera  have  enhanced  m^'s.  Stead  and 
Helmberger  (1985)  argue  that  this  effect  may  be  due  to  scattered 
Rayleigh  waves.  Blandford  and  his  colleagues  have  suggested  this 
complexity  earlier  with  respect  to  in-coming  signals  and  recipro¬ 
city  arguments.  We  will  briefly  discuss  a  rather  simplified 
model  of  the  phenomena  for  Yucca  Flats. 

Basin  geometry  can  have  several  effects  on  seismic  wave 
propagation.  For  example,  the  angle  of  basin  termination  can 
have  a  strong  attenuating  effect  on  surface  waves  crossing  this 
boundary.  At  Yucca  Flats  the  basin  terminates  at  relatively 
shallow  dips,  probably  less  than  30°  ,  which  causes  energy  to 

leak  out  the  bottom.  Figure  1  shows  the  observed  strong  motions 
for  the  FLASK  event  as  recorded  on  the  vertical  component.  All 
three  components  are  available  and  particle  motion  plots  indicate 
that  the  larger  period  secondary  arrivals  are  Rayleigh  waves. 
Note  their  relative  reduction  in  amplitude  at  the  western  hard 
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Figure  1.  Vertical  components  of  the  strong-motion  records  at  roughly 
3  km  with  peak  amplitudes  indicated  above  each  trace.  Note  the  strong 
reduction  in  amplitudes  towards  the  west. 


a 


r 


Figure  2.  Diagram  indicating  the  enclosed  two-dimensional  surface  and 
corresponding  integral  representation. 
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rock  sites.  The  important  problem,  then,  is  that  assuming 
the  surface  wave  energy  is  generated  isotropically  from  the 
source,  where  does  the  energy  go  and  does  it  leave  as  P  and/or  S 
waves?  Part  of  it  may  appear  in  teleseismic  body  wave  coda  as 
suggested  earlier.  There  are  two  important  observations  which 
can  be  noted.  The  first  is  that  there  are  large  complicating 
anivals  following  the  P  and  pP  phases  for  most  Yucca  events. 
The  second  is  that  for  those  events  showing  the  complications , 
the  later  phases  vary  greatly  as  a  function  of  source  position. 

Approach 

In  order  to  model  the  teleseismic  records  using  a  complex 
source  region,  we  derive  a  representation  theorem  similar  to 
Kirchhoff  integration,  see  Scott  fuo  Heimberger  (1985)  We  will 
begin  with  a  simple  2-D  Kirchhoff  method  outlined  in  Figure  2, 
where  KQ  is  the  modified  TK'sscl  function.  The  summation  proce¬ 
dure  is  similar  to  that  discussed  earlier  except  that  we  need 
only  sum  over  a  line  of  elements. 

An  application  of  this  technique  to  an  idealized  cross-sec¬ 
tion,  taken  from  NTS  is  given  in  Figures  3  and  4.  Note  that  the 
receivers  are  arranged  from  roughly  the  location  of  the  Yucca 
Fault  across  the  Valley  and  onto  hard  rock.  The  simplified 
cross-section  of  Yucca  Flat  produces  the  surface  waves  shown  in 
the  accelerograms  of  Figure  3.  As  above,  these  are  computed 
using  a  two-dimensional  finite  difference  method  and  include  a 
Von  Seggern-Blandf ord  RDP  source  with  a  K  of  12  and  B  of  1. 
These  accelerograms  were  produced  using  the  asymmetric  source  of 
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Figure  3.  Upper  portion  displays  a  two-dimensional  cross-section  along 
with  the  locations  of  shots  and  receivers.  The  lower  section  shows  the 
comparison  of  synthetics  (velocity)  for  three  cases. 
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Vidale  and  Helmberger  (1986),  discussed  in  the  next  section. 
The  receivers  remain  fixed  for  all  three  cases,  and  the  eastern¬ 
most  station  in  each  case  is  the  first  receiver  immediately  to 
the  west  of  the  source.  The  first  case  is  a  flat  layer  geometry 
with  the  same  velocities  as  the  simplified  cross-section  and  is 
included  for  comparison.  Comparing  Source  1  and  the  flat  layer 
geometry  shows  an  increase  in  surface  wave  frequency  and  greater 
amplitude  decrease  as  the  waves  pass  into  hard  rock  across  the 
boundary.  Source  4  demonstrates  well  the  drop  in  surface  wave 
amplitude  when  crossing  this  boundary  as  is  observed  for  many 
Yucca  events. 


Teleseismic  Synthetics  -  WWSSN 


Line  explosion  with  correction  term 


a— 20° 


Line  explosion  uncorrected 


Figure  4.  Comparison  of  teleseismic  synthetics  obtained  from 
propagating  the  signals  out  the  bottom  of  the  model  via  Kirchoff.  Note 
the  strong  secondary  arrival  associated  with  case  4. 
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The  teleseismic  synthetics  of  Figure  4  are  computed  using 
the  two-dimensional  Kirchhoff  method.  Once  again  the  flat  layer 
case  is  included  for  comparison.  The  record  sections  show  the 
variation  in  the  waveforms  as  the  source  is  moved  across  the 
basin.  The  most  interesting  observation  is  that  moving  the 
source  one  kilometer  in  the  basin  is  much  more  important  than 
varying  the  take-off  angle  five  degrees.  Also  shown  is  a 
comparison  with  a  symmetric  explosion  demonstrating  the  impor¬ 
tance  of  the  correction  mentioned  above. 

Correction  Factor  for  Explosions  in  2-D  Numerical  Grids 

We  will  show  that  the  amount  of  energy  leaving  the  source 
region  at  an  angle  0  with  the  vertical  in  the  2-D  grid  may  be 
approximated  by  the  amount  of  energy  in  the  point  source  case 
multiplied  by  ^sin  0.  The  additional  sin  0  in  the  point  source 
or  three-dimensional  solution  can  be  explained  in  terms  of 
geometrical  spreading  as  is  shown  in  Figure  5a  and  5b.  The 
energy  with  takeoff  angle  i  between  ig  and  ig  +  dig  for  the  point 
source  becomes 

Ep  a  12a  sin  i _ £_) _ £  diQ 

2  n  r  2 


or 


Ep  a  sin  ig  dig 

while  for  the  line-source  EL  a  dig 

Since  energy  is  proportional  to  the  square  of  the  amplitude  we 

v — ; - 

obtain  the  sin  0  dependence. 

If  we  use  an  isotropic  explosion  as  the  source  in  the  2-D 
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Figure  5.  Diagrams  showing  energy  with  takeoff  angle  i  in  the  range 
i0  <  i  <  i0  +  di  for  a  point  source  and  for  a  line  source.  The  energy 
varys  as  sin  for  the  point  source  case  but  does  not  vary  as  a  function 
of  0  for  the  line  source  case. 


Figure  6.  Comparison  of  radiation  pattern  for  corrected  and  uncorrected 
line  sources.  The  level  line  shows  the  isotropic  radiaion  pattern 
which  results  from  an  uncorrected  explosive  line  source.  The  sin0 
curve  shows  the  best  radiation  pattern  to  simulate  an  explosive  point 
source.  The  two  sinusoidal  curves  show  the  result  of  mixing  a  line 
source  force  with  a  line  source  explosion  with  50%/50%  and  60%/40% 
weighting.  The  mixed  sources  are  meant  to  be  accurate  in  the  range 
0  =  +20°  to  +1600. 
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grid,  each  arrival  in  a  record  may  have  a  different  take-off 


angle  0,  but  we  can  only  correct  for  a  constant  V  sin  6.  The 
result  is  that  the  importance  of  vertically  as  opposed  to 
horizontally  traveling  energy  is  overemphasized  in  the  line 
source  compared  to  the  point  source  case. 

One  might  ask  why  not  simple  multiply  the  isotropic  source 
by  V  sin  0.  Such  a  source  does  not  satisfy  the  2-D  elastic  wave 
equation,  unfortunately,  and  will  not  maintain  the  V  sin  0 
radiation  pattern  once  the  energy  leaves  the  source  region, 
primarily  due  to  the  cusp  in  the  V  sin  0  at  0=0°. 

The  source  functions  we  have  found  which  are  solutions  to 
the  2-D  elastic  wave  equation  have  radiation  patterns  of  sinn  0 
cosm  0,  where  we  may  choose  n  and  m.  An  isotropic  line  source 
explosion,  for  example,  is  the  solution  with  n  =  m  =0,  and  the 
dislocation  sources  have  n  +  m  =  2  (Vidale  and  Helmberger,  1986). 
Also,  because  of  the  asymptotic  nature  of  our  solutions,  the 
compressional  and  shear  parts  of  the  source  separate.  The 
correction  we  introduce  is  to  add  the  compressional  component  of 
the  horizontally-directed  force  term  (n  =  1  and  m  =  0)  to  the 
explosive  source  that  decreases  the  amplitude  of  energy  leaving 
the  source  vertically,  but  leaves  unchanged  the  amplitude  of 
energy  going  horizontally  out  the  side. 

By  judiciously  mixing  the  explosive  and  force  terms,  we  can 
modify  the  vertical  radiation  pattern  of  the  explosion  to  mimic 
the  V  sin  0  we  desire.  Figure  6  shows  the  radiation  patterns  that 
result  from  using  100%  explosion,  50%  explosion  mixed  with  50% 
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force,  and  40%  explosion  mixed  with  60%  force.  These  cases  are 


compared  with  the  Vsin  S'  and  isotropic  line  source  radiation 
patterns.  Energy  that  leaves  the  source  at  angles  near  0  =  90° 
is  not  affected  by  the  correction,  but  energy  at  angles  near  0  = 
0°  is  markedly  affected.  The  mix  of  explosive  and  force  source 
expressions  determines  where  in  the  radiation  pattern  the  source 
is  most  accurate.  The  50%/50%  mix  is  most  accurate  near  0  =  90° 
and  the  60%/40%  mix  is  less  accurate  near  0  =  90°,  but  is  also 
accurate  near  0  =  30°,  as  may  be  seen  in  Figure  6.  It  is  clear 
from  Figure  6  that  only  energy  leaving  the  source  at  positive 
angles  may  be  modeled  with  this  corrected  source. 

Recommendations 

Line-source  calculations  to  generate  point-source  spreading 
with  proper  radiation  patterns  looks  quite  promising.  However, 
the  problem  is  truly  3 -dimensional  and  we  plan  to  perform  some 
full  scale  checks  in  the  near-future.  The  final  objective,  of 
course,  is  relating  the  above  synthetics  to  actual  data  for  known 
geology  and  the  proper  geometries. 
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MODELING  Lg  PROPAGATION  ACROSS  CRUST-MANTLE  TRANSITION  ZONES 
AND  EVALUATING  SURFACE  WAVE  PATH  CORRECTIONS 


Janice  Regan,  Peter  Glover,  and  David  G.  Harkrider 

Seismological  Laboratory,  California  Institute  of  Technology, 

Pasadena,  California,  91125 

INTRODUCTION 

Finite  element  forms  of  the  seismic  Representation  Theorem  (RT)  have  been  developed  to  allow 
modeling  of  teleseismic  Lg  waves  crossing  continental  margins,  or  other  heterogeneities.  The  accuracy  of 
this  formulation  and  the  efficiency  of  an  absorbing  boundary  condition  designed  to  reduce  computation 
time  have  been  previously  established.  Initial  calculations  of  Lg  wave  propagation  through  transition 
regions  have  been  verified  and  extended  by  a  numerical  experiment  which  studied  the  effects  of  different 
lengths  of  transition  regions  on  the  Lg  waves  traversing  them.  The  FE  calculations  are  driven  by  specify¬ 
ing  the  displacement  time  histories  at  a  column  of  nodes  which  represent  the  boundary  between  the 
laterally  homogeneous  region  and  the  transition  zone.  These  time  histories  can  be  calculated  using  stan¬ 
dard  techniques  when  the  source  is  separated  from  the  margin  by  long  continental  paths.  Alternatively, 
they  can  be  input  from  other  FE  calculations.  Both  of  these  approaches  are  used  in  the  results  reported 
here.  The  first  approach,  will  be  referred  to  as  the  hybrid  method.  The  time  histories  used  as  FE  input 
for  the  hybrid  calculations  described  below  are  Lg  mode  sum  displacement  time  histories  calculated  at  .5 
km  vertical  spacing  down  a  grid  edge  at  a  distance  a=  1500  km  from  a  8  km  deep  source.  The  FE  cal¬ 
culation  propagates  the  energy  through  the  remaining  portion  of  the  path,  a  distance  P,  to  the  receivers 
located  at  distance  6.  We  should  point  out  however  that  our  technique  is  only  valid  for  waves  propagating 
at  normal  incidence. 

Recent  estimates  of  explosion  moments  using  surface  waves  are  made  using  surface  wave  path 
corrections  based  on  a  modification  of  a  two  dimensional  approximation  of  normal  incident  surface  waves 
crossing  slowly  varying  transition  zones.  Similar  to  this  correction,  the  Gaussian  beam  approximation  for 
surface  waves  also  conserves  lateral  energy  flux  .  In  addition  it  includes  the  effect  of  lateral  heterogenity 
on  geometrical  or  optical  spreading.  Since  neither  technique  was  intended  for  generating  path  corrections 
for  surface  waves  generated  by  sources  in  limited  source  regions  with  possibly  sharp  boundaries,  such  as 
low  velocity  basins,  it  is  important  to  test  their  validity.  This  was  done  by  comparing  synthetic  Rayleigh 
waves  calulated  by  these  techniques  with  those  obtained  from  hybrid  finite-element  codes  coupled  to  ana¬ 
lytic  teleseismic  modal  propagator  codes. 

Lg  ACCOMPLISHMENTS 
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An  initial  test  of  the  effects  of  a  transition  region  on  the  propagation  of  Lg  mode  sum  seismograms 
used  a  pair  of  20  km  transitions  of  the  type  shown  in  Figure  1.  The  results  of  this  test,  outlined  below, 
indicated  that  a  numerical  experiment  to  quantify  the  dependence  of  Lg  attenuation  across  transition 
regions  on  the  structure  of  those  transition  regions,  could  give  some  significant  insights  into  where  and 
why  Lg  disappears  when  it  interacts  with  oceanic  structures.  The  initial  test  results  also  describe  the 
common  properties  observed  for  all  transitions  studied.  Energy  reflected  back  from  the  transition  zone 
toward  the  source  is  clearly  visible.  The  energy  impinging  on  the  crust  to  water  sloping  boundary  is 
mostly  reflected.  The  passage  of  the  Lg  wave  through  the  forward  transition  region  effected  the  distribu¬ 
tion  of  energy  with  depth.  Since  no  energy  is  transmitted  into  the  water  layer,  the  remaining  energy  must 
be  concentrated  by  the  sloping  boundary  into  the  oceanic  crustal  layer.  Thus,  the  amplitudes  in  the  oce¬ 
anic  crustal  layer  at  the  end  of  the  transition  are  larger  than  the  amplitudes  between  5  km  and  ten  km 
depth  at  the  beginning  of  the  transition.  The  difference  is  maximum  at  a  depth  of  5km  and  decreases 
with  depth,  suggesting  that  some  energy  is  traveling  down  the  sloping  surface  into  the  oceanic  crustal 
layer.  The  energy  incident  on  the  sloping  crust  to  mantle  boundary  seems  to  be  primarily  transmitted. 
Conversely,  for  the  reverse  transition,  the  energy  traveling  in  the  oceanic  crustal  layer  spreads  into  the 
entire  continental  crustal  layer  reducing  the  amplitudes  at  any  particular  depth.  Most  of  the  energy 
incident  on  the  sloping  mantle  to  crust  boundary  is  reflected.  As  the  wavefront  goes  through  the  con¬ 
tinent  to  ocean  transition  surface  displacement  history  amplitudes  increase  as  the  crustal  thickness 
decreases.  When  the  wavetrain  propagates  through  the  oceanic  structure  surface  amplitudes  decay  with 
distance,  rapidly  at  first,  and  then  more  slowly.  Furthermore,  as  the  wavefront  passes  through  the  oce¬ 
anic  to  continental  transition  a  rapid  decay  in  amplitude  is  observed.  After  the  Lg  wave  reenters  the  con¬ 
tinental  structure  amplitude  decay  continues  for  a  short  distance  while  the  waveform  and  amplitude  sta- 
blize. 

A  study  of  the  effects  of  the  length  of  simple  transition  regions  on  the  attenuation  of  SH  type  Lg 
mode  sum  seismograms  passing  through  them  has  yielded  some  interesting  results.  In  order  to  discuss 
these  results  the  models  must  first  be  described  and  the  methods  of  analysis  explained.  Two  classes  of 
transition  models  were  considered.  An  example  of  each  class  is  illustrated  in  Figure  1.  Calculations  were 
performed  for  four  individual  models  from  each  class,  for  a  continental  crustal  layer  over  a  mantle  half 
space,  and  for  an  oceanic  crustal  layer  over  a  mantle  half  space.  The  difference  between  individual  transi¬ 
tion  models  was  the  length  of  the  transition  region,  or  the  horizontal  distance  between  points  B  and  D 
shown  in  Figure  1  As  the  length  of  the  transition  increases  the  angle  that  the  the  ocean  to  crust  boun¬ 
dary  or  crust  to  ocean  boundary  makes  with  the  horizontal  (<p0(.  or  4>cg  in  Figure  1)  varies  between  3“ 
and  90“  ,  and  the  angle  the  the  crust  to  mantle  boundary  or  the  mantle  to  crust  boundary  makes  with 
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the  horizontal  (<fime  or  4>cm  in  Figure  1)  varies  from  12°  to  90°  .  The  differences  in  slope  of  the  boundaries 
and  the  different  elastic  properties  of  the  layers  the  separate  indicate  that  different  behavior  should  be 
expected  along  those  two  boundaries.  The  initial  model  discussed  previously  used  a  transition  length  of 
20  km.  Real  ocean  to  continent  type  transitions  occur  over  lengths  of  order  100  km.  Limiting  FE  compu¬ 
tation  time  per  model  to  approximately  one  cpu  day  imposed  an  upper  limit  of  100  km  on  the  length  of 
the  transition.  Thus,  the  lengths  used  for  this  investigation  were  a  step  transition  (0  km),  25  km,  50  km. 
and  100  km.  For  each  model  seismograms  were  recorded  at  intervals  of  approximately  5  km  along  the 
surface.  Figure  2  shows  the  RMS  surface  amplitudes  calculated  over  the  first  60  seconds  of  seismogram 
for  each  of  the  receivers  along  the  surface  of  the  50  km  transition.  Seismograms  were  also  recorded,  with 
vertical  spacing  of  2.5  km,  at  several  distances.  Depth  sections  were  recorded  at  distances  A  through  E 
(Figure  l)  for  all  models.  For  the  100  km  transitions  additional  sections  were  recorded  midway  between  B 
and  C  and  midway  between  C  and  D.  For  the  50  km  forward  transition  an  additional  depth  section  with 
vertical  spacing  of  .5  km  was  recorded  25km  beyond  the  end  of  the  transition.  This  section  was  used  as 
input  for  later  calculations. 

The  first  class  of  models  are  of  continent  to  ocean,  or  ’forward’  transitions.  Each  ’forward’  transi¬ 
tion  model  uses  the  same  set  of  input  data,  or  forcing  functions.  In  each  case  the  leftmost  column  of  nodes 
of  the  forward  transition  FE  grid  are  constrained  to  move  with  the  displacement  time  histories  specified 
by  the  driving  functions.  In  each  forward  transition  model  the  transition  region  is  characterized  by  a  con¬ 
tinuous  rate  of  thinning  of  the  crustal  layer  between  the  thirty  two  km  thick  continental  crust  at  the 
beginning  of  the  transition  region  and  the  5  km  oceanic  crust,  overlain  by  5  km  of  water,  at  the  end  of 
the  transition  region.  The  mode  sum  forcing  functions  are  also  used  as  input  to  a  reference  model  of  a 
thirty  two  km  continental  layer  over  a  mantle  half  space.  Depth  sections  at  A  and  B  (see  Figure  la)  were 
subtracted  from  depth  sections  recorded  after  the  same  distance  of  FE  propagation  in  the  continental 
reference  model.  These  differences  give  the  wave  reflected  back  from  the  transition  region. 

The  second  class  of  models  are  of  ocean  to  continent,  or  ’reverse’  transitions.  The  driving  functions 
for  the  ’reverse’  transition  tests  are  recorded  during  the  50  km  ’forward’  transition  calculation.  They  con¬ 
sist  of  a  depth  section  of  hybrid  seismograms  recorded  25  km  past  the  oceanic  end  of  the  50  km  forward 
transition,  which  corresponds  to  a  distance  of  1755  km  from  the  source.  The  vertical  spacing  within  the 
depth  section  is  .5  km.  The  reverse  transition  driving  functions  are  also  used  as  input  to  an  oceanic  refer¬ 
ence  model  of  a  5  km  crustal  layer  over  a  mantle  half  space.  Each  reverse  transition  is  modeled  as  a 
smooth  increase  in  thickness  of  the  crustal  layer  between  a  5  km  thick  oceanic  crust,  overlain  by  5  km  of 
ocean,  at  the  beginning  of  the  transition  and  a  32  km  continental  crust  at  the  end  of  the  transition. 
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The  ’forward’  transition  tests  illustrate  several  properties  that  change  with  the  length  or  slope  of  the 
transition  model.  At  the  crust  to  water  boundary  the  amplitude  of  the  waveform  reflected  back  into  the 
continental  layer  ,  that  is  back  toward  the  source,  decreases  as  the  length  of  the  transition  increases.  A 
similar  effect  of  much  smaller  magnitude  is  observed  along  the  crust  to  mantle  boundary.  The  difference 
in  the  magnitudes  of  the  effects  at  these  two  boundaries  can  be  explained  by  the  difference  in  velocity 
contrast  across  them  and  the  difference  in  the  angles  4>mc  and  <j)c0  .  Although  no  energy  in  transmitted 
into  the  water  layer,  not  all  energy  impinging  on  the  water  layer  is  reflected.  Instead  some  of  the  energy  is 
concentrated  by  the  transition  into  the  oceanic  crustal  layer.  The  amount  of  concentrated  energy, 
increases  as  the  length  of  the  transition  increases  (Figure  3)  .  Concentrated  energy  is  observed  as  an 
increase  in  amplitude  of  seismograms  in  the  oceanic  crustal  layer  with  respect  to  the  amplitudes  of  seismo¬ 
grams  at  the  same  depth,  between  5  km  and  10  km,  in  the  continental  crustal  layer  at  the  beginning  of 
the  transition.  Due  to  the  transmission  of  energy  across  the  crust  to  mantle  boundary,  detection  of  a  simi¬ 
lar  concentrating  effect  associated  with  this  boundary  will  require  more  detailed  analysis  of  the  calculated 
results.  Amplitudes  in  the  mar'  e  layer  within  the  transition  region  are  maximum  immediately  below  the 
crust  to  mantle  transit--)”  „nd  decay  nearly  exponentially  with  depth  below  the  boundary.  The  ampli¬ 
tude  of  the  transmitted  seismograms  below  the  oceanic  crustal  layer,  at  the  oceanic  end  of  the  transition, 
decrease  as  th.  iength  of  the  transition  increases.  Figure  2  shows  that  the  seismogram  amplitudes  at  the 
surface  of  the  crustal  layer  show  a  general  increase  as  one  passes  through  the  transition,  however,  this 
increase  does  not  appear  to  be  monotonic.  The  variation  in  surface  amplitude  within  the  transition  zone 
has  similar  pattern  for  all  lengths  of  transitions,  but  the  width  of  the  oscillations  scale  as  the  width  of  the 
transition  region.  This  is  illustrated  in  Figure  3,  which  shows  a  sketch  of  the  variation  in  the  behavior  of 
RMS  amplitudes  with  distance  along  the  surface  for  each  transition  length. 

The  ’reverse’  transition  tests  have  not  been  as  extensively  studied.  However,  it  is  clear  that  a  longer 
transition  is  much  more  effective  at  reducing  surface  amplitudes.  Examination  of  Figures  2  and  3  show 
that  the  surface  amplitude  decay  across  the  transition  is  monotonic  for  the  reverse  transition,  and  that  the 
surface  amplitude  vs.  distance  relation  along  the  oceanic  portion  of  the  path  is  oscillatory.  This  oscilla¬ 
tory  behavior  is  also  observed  for  the  oceanic  reference  model  calculation,  where  it  is  superimposed  upon 
a  slow  decrease  in  amplitude.  The  amplitude  distributions  with  depth  indicate  that  energy  is  being  chan¬ 
neled  down  into  the  mantle  by  the  reverse  transition.  The  amount  channeled  energy  seems  to  increase  as 
the  length  of  the  transition  increases.  The  amplitudes  below  the  crustal  layer  of  the  transition  region  do 
not  die  off  as  strongly  with  depth  as  they  do  in  the  ’forward’  transition  case.  It  appears  that  it  is  in  the 
’reverse’  transition  and  the  propagation  through  the  intermediate  oceanic  layer  that  causes  the  attenua¬ 
tion  of  Lj  type  waves  when  they  travel  across  a  mixed  path  containing  an  oceanic  region. 
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The  length  of  the  intermediate  oceanic  path  between  the  continent  to  ocean  and  the  ocean  to  con¬ 
tinent  transition  is  also  being  investigated.  Calculations  for  three  distances  have  been  done  using  multiple 
FE  calculations  to  go  from  the  edge  of  the  first  transition  model  to  the  receiver.  Due  to  computational 
limits,  the  lengths  of  these  intermediate  paths  are  short,  but  definitely  within  the  range  of  interest. 
Lengths  of  30,  70,  and  120  km,  have  been  used.  For  the  shortest  paths  energy  transmitted  in  to  the  man¬ 
tle  layer  has  little  time  to  dissipate  out  of  the  calculation  before  much  of  it  passes  back  into  the  crustal 
layer  as  it  travels  through  the  reverse  transition.  This  implies  that  there  may  be  a  critical  distance  where 
enough  of  the  energy  has  been  lost  by  propagation  through  the  oceanic  structure  that  Lg  would  essen¬ 
tially  disappear  due  to  dilution  caused  by  propagation  through  the  reverse  transition  region.  Further 
analysis  of  completed  calculations,  and  additional  calculations  are  needed  to  establish  whether  this  is  true, 
and  to  quantify  any  such  effect.  In  particular  a  coupling  method  to  pass  the  FE  results  back  into  the 
mode  sum  calculations  is  nearing  completion.  This  will  allow  consideration  of  longer  paths,  and  easier 
determination  of  shorter  paths. 

SURFACE  WAVE  PATH  CORRECTION  ACCOMPLISHMENTS 

The  source  region  medium  is  a  cylinder  with  a  radius  of  1.8  km  and  extends  from  the  free  surface 
down  to  a  depth  of  1.8  km.  The  source  cylinder  is  embedded  in  a  half  space  of  the  same  elastic  properties 
as  that  of  the  upper  part  of  the  propagating  medium.  The  RT  surface,  over  which  the  propagation 
Green’s  functions  are  integrated,  surrounded  the  plug  at  a  distance  of  0.3  km  into  the  half  space  of  propa¬ 
gation  material.  The  RT  surface  is  composed  of  11  nodes  or  rings  per  side  with  a  spacing  of  0.2  km.  The 
propagation  model  is  CIT109,  which  is  based  on  western  US  observations  and  has  a  14  km  thick  upper 
crust.  The  three  elastic  materials  used  in  the  source  medium  cylinder  are  characteristic  of  the  shot  point 
properties  at  three  NTS  areas;  Climax  Stock,  Pahute  Mesa  and  Yucca  Flat.  The  seismic  velocities,  densi¬ 
ties  and  other  material  properties  for  the  source  media  are  given  in  the  following  table. 


Material  Properties  of  Source  Region 

Models 

SOURCE 

MODEL  (S) 

a 

km/s 

0 

km/s 

VCIT 

Vs 

(X+2  h)cit 
(\+2n)s 

a2/? 

YUCCA 

2.35 

1.30 

1.86 

10.7 

10.2 

3.27 

PAHUTE 

4.00 

1.90 

2.30 

4.06 

2.86 

4.43 

CLIMAX 

5.33 

2.78 

2.67 

1.63 

1.39 

3.68 

CIT109 

6.20 

3.51 

2.74 

1.00 

1.00 

3.12 

In  order  to  demonstrate  that  modeling  an  explosion  by  applying  a  stress  glut  at  two  elements  along 
the  the  vertical  axis  was  sufficient,  we  made  an  RT  calculation  in  which  the  source  medium  was  the  same 
as  the  surrounding  medium.  We  then  calculated  its  analytic  Rayleigh  wave  analog,  a  point  explosion  at  a 
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depth  of  0.4  km  in  the  CIT109  model.  The  time  domain  records  of  both  calculations  as  seen  through  a 
long  period  VVWSSN  instrument  at  3000  km  were  essentially  indistinguishable.  The  maximum  error  was 
less  than  2  %. 

Previous  reports  have  focused  on  time  domain  comparisons  between  RT  and  CLEF.  Since  the  phase 
of  the  CLEF  calculation  is  not  used  to  estimate  the  explosion  and  tectonic  release  moments,  we  will  res¬ 
trict  this  report  to  comparisons  of  amplitude  spectra  and  spectral  ratios.  In  Figure  4a,  we  show  the  spec¬ 
tra  and  spectral  ratios  of  the  analytic  and  RT  synthetics  for  the  pure  CIT109  model  with  the  WWSSN-LP 
instrument.  The  calculation  frequency  range  is  0  to  2.5  hz.  The  agreement  between  the  two  calculations  is 
excellent  out  to  the  roll-off  frequency  of  the  cosine  filter  which  was  applied  to  the  finite  element  forcing 
functions  before  their  decimation  to  0.2  second  time  steps.  The  original  finite  element  time  step  was  con¬ 
siderably  smaller,  by  nearly  a  factor  of  13. 

In  this  figure,  we  only  spectra  out  to  0.2  hz.  This  frequency  range  gives  the  best  detail  of  the  fre¬ 
quencies  of  interest  in  estimating  seismic  moments  using  surface  wave  path  corrections.  At  periods  greater 
than  50  seconds,  the  RT  calculation  has  larger  amplitudes  than  the  analytic  or  exact  synthetic.  This  can 
be  considered  the  long  period  upper  limit  to  the  accuracy  of  the  RT  calculations. 

In  the  rest  of  Figure  4,  we  display  the  spectra  and  spectral  ratios  of  RT  calculations  for  the  three 
source  media  and  their  CLEF  approximations.  The  agreement  between  RT  and  CLEF  spectra  is  good 
especially  considering  the  rational  behind  the  approximations.  The  three  cases  all  show  the  same  charac¬ 
ter  in  that  CLEF  spectra  is  smaller  than  RT  at  long  periods  and  greater  at  higher  frequencies.  The  cross¬ 
over  is  somewhere  between  20  and  10  seconds  period. 

In  Figure  5  a  and  h,  we  show  the  spectra  of  the  RT  model  of  Yucca  with  a  latteral  and  vertical  tran¬ 
sition  zone  out  to  0.2  hz  and  1.0  hz  respectively.  It  is  compared  with  its  CLEF  analogue.  Beyond  0.2  hz, 
the  CLEF  approximation  becomes  increasingly  bad.  This  is  true  of  all  the  models  with  the  CLEF  approxi¬ 
mation  having  higher  frequency  content.  Because  of  the  small  vertical  extent  of  the  inhomogeneous 
source  region  relative  to  the  wavelengths  of  the  surface  waves,  the  difference  between  the  spectra  of  CLEF 
and  GBSW  is  negligible  for  periods  greater  than  two  to  four  seconds.  This  is  shown  in  Figure  5  c  for  the 
Yucca  models.  At  higher  frequencies,  the  GBSW  approximation  has  less  frequency  content  out  to  1.0  hz  as 
shown  in  Figure  5  d. 

In  Figures  6  and  7,  we  show  the  CLEF  and  GBSW  approximations  out  to  1.0  hz  compared  to  the 
RT  models  discussed  above.  The  CLEF  comparisons  are  on  the  left  and  the  GBSW  comparisons  are  on 
the  right  of  each  figure.  In  every  model  calculated,  the  GBSW  outperforms  the  CLEF  approximation  and 
never  has  less  frequency  content  than  the  RT  model. 


32 


rwari  Tnfcftitkton 


'*  l hr  iT4n>iiittn  rcfitott.  fraclivit'-  uf  t hr  innsitiuh  Irngfii  All  j. 

Ai'rj-i  fur  Miiitt  ib  ilit  100  km  (r*iiiiiioo  rctulu 


GEOPHYSICAL  INVESTIGATIONS  AT  PAHUTE  MESA .  NV 


JOHN  F.  FERGUSON 

THE  UNIVERSITY  OF  TEXAS  AT  DALLAS 
This  paper  is  a  progress  report,  midway  though  a  two  year 
study  a-f  the  shallow  (<5  km)  geophysical  structure  under  Pahute 
Mesa,  NV  and  its  e-f-fect  on  teleseismic  observations  of  events  at 
that  site.  Topics  to  be  discussed  include  the  geophysical 
properties  as  deduced  from  well  logs,  new  gravity  stations 
accquired  in  June  1985,  models  of  the  Bouguer  gravity  for  the 
Silent  Canyon  Caldera  and  a  seismic  reflection  -  refraction 
experiment  to  be  conducted  in  May  and  June  of  1986. 

Both  uphol e  P-wave  velocity  surveys  and  either  density  or 
borehole  gravity  logs  were  available  in  holes  U19af ,  Uel9z,  U19x, 
U19ak,  U19ae  and  U19ab.  Locations  for  these  holes  are  shown  in 
Figure  1.  The  geologic  formations  can  be  geophysically  clustered 
into  four  units  as  follows.  1)  Tuffs  -  low  density  and  velocity. 
2)  Lavas  -  high  density  and  velocity.  3)  Rainer  Mesa  tuff  (Tmr)  - 
high  density  and  law  velocity.  4)  Peralkaline  rocks  of  the  Silent 
Canyon  Caldera  -  high  density  and  velocity,  which  constitutes  the 
geophysical  "basement".  The  tuffs  have  been  characteri zed  by 
linear  density  and  velocity  functions  as  illustrated  in  Figure  2. 
Below  the  water  table,  at  about  700  m,  the  gradient  decreases  to 
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0.42  g m/cm  /km,  with  similar  behavior  in  velocity  also  expected. 

Prior  to  1985  the  gravity  data  base  for  the  Pahute  Mesa  - 
Timber  Mt,  area  consisted  of  3200  stations  obtained  by  the  USGS 
in  the  early  1963's.  Elevation  control  was  based  on  poor  maps  and 
barometric  altitudes.  The  station  density  was  generally  1  or  2 
/km*".  In  june  of  1985  UTD  obtained  300  new  stations  an  traverses 
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The  elevation  control  was  by  EDM, 


along  reads  at  300  m  intervals, 
so  that  the  new  data  is  approximately  an  order  o-f  magnitude  more 
accurate  than  the  old.  Near  the  southeast  margin  o-f  the  Silent 
Canyon  Caldera,  two-dimensional  coverage  was  obtained.  In  Figure 
3  the  19B5  traverses  are  shaded  over  the  station  locations.  The 
old  and  new  data  agree  well  in  areas  o-f  overlap.  The  Bouguer 
anomaly  was  computed  far  the  merged  data  with  a  density  of  2.0 
gm/cm 3  and  terrain  corrections  were  applied  using  digital 
topographic  data.  The  resulting  Bouguer  anomaly  is  contoured  in 
Figure  4. 

The  25  mGal  residual  anomaly  can  be  modeled  as  a  low  density 
(with  linear  density  function)  surface  inclusion  of  volcanic  fill 
in  the  high  density  Silent  Canyon  collapse  structure.  The  high 
density  Tmr  and  lava  units  are  well  controlled  by  boreholes  and 
their  effect  can  be  stripped  from  the  data.  The  Silent  Canyon 
interface  is  more  sparsely  controlled,  but  can  be  estimated  by 
inversion  using  a  modification  of  Oldenburg's  (1974)  method.  The 
regional  gravity  anomaly  is  simultaneously  estimated  by  a  least 
squares  fit  to  borehole  constraints.  As  an  example  the  profile  B- 
B'  on  Figure  4  is  presented  in  Figure  5.  The  Tmr  is  shown  by  the 
dashed  pattern  and  the  lava  unit  (Tra)  is  presented  with  the 
cuneform  pattern.  The  Silent  Canyon  basement  is  stippled. 
Excellent  agreement  is  found  at  four  boreholes  near  the  profile. 
The  estimated  regional  shows  a  small  "high",  which  may  correspond 
to  the  deep  crustal  high  velocity  anomally  noted  by  previous 
investigators  (Taylor,  19S3)  . 

In  May  and  June  of  1986  a  seismic  reflection  -  refraction 
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experiment  will  be  performed  at  Pahute  Mesa  by  UTD  in  cooperation 
with  Los  Alamos  National  Laboratory  and  the  USGS  re-fraction  group 
at  Menlo  Park,  CA.  Seventeen  high  explosive  shots  will  be  fired 
in  outcropping  Silent  Canyon  basement  at  the  location  indicated 
in  Figure  4  and  recorded  on  the  line  crossing  the  southeast 
caldera  margin  and  trending  northwest  for  15  km.  Up  to  70 
Dinoseis  shots  will  be  recorded  at  near  offsets.  A  24  channel, 
floating  point  recording  system  will  be  used  with  60  m  geaphone 
group  interval  and  1400  m  spreads  with  50%  overlap. 

The  USGS  (Hoffman  and  Mooney,  1984)  has  had  considerable 
success  with  a  combined  refraction  and  gravity  i ntepretat i on  at 
the  Crater  Flat  Caldera,  to  the  south  of  Pahute  Mesa.  The  model 
presented  in  Figure  6  was  derived  from  data  that  are  more  sparse 
than  the  Pahute  Mesa  experiment,  but  with  much  longer  offsets 
(out  to  40  km)  for  deeper  crustal  penetration.  The  1986 
experiment  will  map  the  Silent  Canyon  basement  and  shallow 
structures  in  the  tuffs  and  lavas  filling  the  caldera.  It  may  be 
possible  to  record  reflections  from  as  deep  as  5  km,  which  could 
be  associated  with  Paleozoic  units  and  Tertiary  or  Mesozoic 
i ntrusi ves. 

In  the  future  it  would  be  passible  to  extend  the  depth 
range  of  the  model  by  recording  nuclear  tests  at  Pahute  Mesa, 
Rainer  Mesa  and  Yucca  Flat  with  event  recorders  on  the  same 
profile.  This  concept  is  illustrated  in  Figure  7.  Detailed 
knowledge  of  the  shallow  structure  at  Pahute  Mesa  gained  in  1986 
and  previously  far  Rainer  Mesa  and  Yucca  Flat  will  permit 
i n terpr etat i on  of  the  far  offset  data  for  mid— crustal  structure. 

I  would  like  to  acknowledge  the  contributions  of  students  at 
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UTD  to  this  project  and  especially  Sharon  Reamer  for  her  part  in 
the  gravity  investigations.  Much  a-f  this  study  would  have  been 
impossible  without  the  support  a-f  personnel  at  Los  Alamos 
National  Laboratory.  Allen  Cogbill,  Rick  Warren,  Jack  House  and 
Don  Collins  a-f  Los  Alamos  deserve  special  recognition. 

REFERENCES 

Ho-f-fman,  L.  R.  ,  (1984),  A  seismic  study  o-f  Yucca  Mountain  and 

vicinity,  southern  Nevada;  data  report  and  preliminary  results, 
USGS  Open-File  Report  83-588,  Menlo  Park,  CA. 

Oldenburg, D.  W. ,  (1974),  The  inversion  and  i n terpretat i on  of 

gravity  anomalies,  Geophys. ,  v  39,  p  526-536. 

Taylor,  S.  R. ,  (1983),  Three-dimensional  crust  and  upper  mantle 

structure  at  the  Nevada  Test  Site,  J.  Geophys.  Res. ,  v  88,  p2220— 


1  ’  8  0  37’  30’ 
3 7  0  2 2'  30’  + 


N 


37°  7’  30’  +■ 


1  16°  15’ 


Figure  t 


39 


Velocity: 

Tuffs- 

v(z)  =  1.70  +  2.15  z 
Lavas- 

v(z)  =  3.50  km/soc 
Density: 

Tuffs  (excluding  Tmr)- 

p(z)  =  1.51  +  0.74  z 
Lavas  and  Tmr- 

p(z)  =  2.20  gm/cm5 


km/sec 

z  in  km 

gm/cm3 


PAHUTE  MESA  GEOPHYSICAL  LOGS 


Figure  2 
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PAHUTE  MESA  GRAVITY  MODEL 


Figure  6 

if  tPr 


Crustal  section  nodlfled  f  roe  Snyder  and  Carr  (1982)  with  layer 
densities  and  velocities  (Table  1),  observed  and  calculated 
travel ties  delays  for  P-wave  arrivals  froa  nuclear  sources,  and 
calculated  and  observed  residual  gravity.  Observed  and  calculated 
delay  tlsii  agree  well,  except  over  Crater  Plat  where  the  observed 
delay  Is  0.2  s  larger  than  that  calculated  (see  text).  Orill  holes 
are  projected  onto  the  profile. 
uoff^an  and  Mooney  (1984) 


after  Taylor  ( 1983) 
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REGIONAL  STUDIES  WITH  BROADBAND  DATA 


T.  V.  McEvilly  and  L.  R.  Johnson 
University  of  California,  Berkeley 

VELOCITY  MODEL  FOR  PAHUTE  MESA 

Considerable  attention  has  been  given  to  the  task  of  estimating  a  velocity 
model  for  the  shallow  structure  at  Pahute  Mesa,  because  such  a  model  is  impor¬ 
tant  for  other  parts  of  the  research  effort.  Data  from  all  of  the  events  that  we 
have  recorded  at  Pahute  Mesa  were  used  in  this  study  (see  Figure  1).  Since  most 
of  these  events  had  depths  near  600  meters,  it  was  possible  to  combine  the 
seismograms  recorded  from  different  events  and  form  a  composite  seismic  section 
for  Pahute  Mesa.  With  the  use  of  this  composite  section,  the  arrival  times  of 
both  P  and  S  waves  were  estimated.  These  arrival  times  were  converted  to  velo¬ 
city  models  by  both  a  trial-and-error  method  and  by  a  formal  inverse  procedure. 
In  both  cases  the  results  were  similar,  and  the  final  model  is  shown  in  Figure  2. 
In  order  to  obtain  information  on  possible  reflecting  horizons  within  the  caldera, 
the  composite  section  was  subjected  to  some  of  the  standard  processing  tech¬ 
niques  that  are  applied  to  controlled-source  reflection  data.  No  definitive 
reflectors  were  revealed  by  this  process. 

MOMENT  TENSOR  INVERSIONS 

The  availability  of  the  refined  velocity  model  for  Pahute  Mesa  allows 
improved  estimates  of  moment  tensors  for  explosions  in  that  region.  The  effort 
has  concentrated  on  two  events  for  which  the  sets  of  near-distance  recordings  are 
most  complete,  Harzer  and  Chancellor.  The  locations  of  these  events  and  their 
recording  stations  are  shown  in  Figure  1  and  the  shot  parameters  are  given  in 
Table  1.  The  station  distances,  azimuths,  and  elevations  are  listed  in  Tables  2 
and  3  along  with  the  maximum  accelerations.  Ground  accelerations  were 
recorded,  but  for  the  purposes  of  the  moment  tensor  inversion  these  were  con¬ 
verted  to  velocity  records,  which  are  shown  in  Figures  3  and  4. 

The  data  in  Figures  3  and  4  were  inverted  to  obtain  estimates  of  the  six  ele¬ 
ments  of  the  symmetric  second-rank  moment  tensor.  The  velocity  model  shown 
in  Figure  2  was  used  to  generate  the  Green  functions.  The  isotropic  part  of  the 
moment  tensor,  which  is  equivalent  to  the  reduced  displacement  potential  of  a 
spherically  symmetric  explosive  source,  is  shown  for  the  two  events  in  Figure  5. 
This  isotropic  part  is  similar  for  the  two  exp  osions,  although  it  is  slightly  smaller 
and  has  a  faster  rise  time  for  the  Chancellor  event.  Both  show  overshoot,  but  the 
result  for  Chancellor  has  a  pronounced  minimum  at  about  2  sec.  This  may  be 
related  to  spall  closure. 

The  deviatoric  part  of  the  moment  tensor  has  also  been  examined.  It  has  a 
lower  frequency  content  than  the  isotropic  part,  seems  to  be  more  variable,  and  is 
less  reliably  determined. 
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SOURCE  LOCATION  IN  A  LATERALLY  HETEROGENEOUS  EARTH 


A  method  has  been  developed  for  locating  seismic  sources  in  a  laterally 
heterogeneous  earth  that  makes  uses  of  empirical  tectonically  regionionalized  P- 
wave  tau  functions  and  slowness-dependent  source  and  receiver  corrections.  The 
method  has  been  tested  by  applying  it  to  28  events  at  the  Nevada  Test  Site 
which  have  known  locations.  Using  only  readings  in  the  interval  15-95  deg,  the 
method  yields  a  mean  mislocation  of  3.32  km  and  a  mean  sample  variance  of 
travel-time  residuals  of  0.32  sec2.  This  was  compared  to  the  standard  method  of 
location  which  assumes  a  laterally  homogeneous  model  and  Jeffreys-Bullen  travel 
times  which  yielded  a  mean  mislocation  of  7.33  km  and  a  mean  sample  variance 
of  0.52  sec2. 

The  method  was  also  applied  to  49  events  in  Eastern  Kazakh  where  the 
mean  sample  variance  was  0.34  sec2  as  compared  to  0.60  sec2  for  the  standard 
method.  For  5  events  in  Novaya  Zemlya  the  mean  sample  variance  was  0.43  sec2 
as  compared  to  0.51  sec2  for  the  standard  method.  For  12  events  in  the  Tuamotu 
Archipelago  the  method  gave  a  mean  sample  variance  of  0.35  sec2  as  compared  to 
0.33  sec2  for  the  standard  method. 
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Table  1. 

Shot  Parameters 

Harzer 

Chancellor 

Date 

6  June  1981 

1  Sept  1983 

Origin  Time 

18h  00m  00.1s 

14h  00m  00.1s 

Latitude 

37.303N 

37.273N 

Longitude 

116.326W 

116.355W 

Surface  elevation 

2097  m 

2040m 

Depth  to  shot 

637m 

625m 

Shot  Medium 

tuff 

tuff 

Depth  to  standing  water 

668m 

647m 

Magnitude  (ML) 

5.4 

5.3 

Table  2.  Network  Data  for  Harzer 

Station 

Distance 

Azimuth 

Elevation 

First  Peak  Acceleration 

Maximum  Acceleration 

(km) 

(deg  E  of  N) 

(m) 

(cm/sec7) 

(cm/sec2) 

Z  R  T 

Z  R  T 

H9 

2.37 

345 

173 

115 

59 

197 

278 

281 

H4 

3.47 

137 

2134 

147 

19 

6 

315 

139 

140 

H2 

3.52 

28 

2060 

121 

115 

-11 

410 

372 

296 

Hi 

4.65 

5 

2112 

123 

59 

-10 

148 

175 

166 

H3 

4.74 

69 

2057 

71 

8 

-2 

163 

257 

144 

H6 

5.61 

206 

2006 

75 

62 

-39 

276 

300 

306 

H8 

5.78 

310 

2060 

37 

4 

-5 

75 

132 

58 

H5 

6.62 

167 

2112 

53 

18 

24 

105 

103 

90 

Table  3.  Network  Data  for  Chancellor 


Test  Site  Boundary 


Figure  1.  Explosion  sites  (open  circles  with  names'}  and  respective  recording  sites  (filled  circles') 
for  near  field  data  acquisition  at  i'abute  Mesa. 
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Figure  2.  Estimates  of  P  and  S  velocities  in  the  shallow  crust  at  Pahute  Mesa. 
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Figure  3.  Velocities  of  ground  motion  at  eight  different  sites  for 
the  event  Harzer.  The  records  have  been  scaled  hy 
epicentral  distance. 
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Figure  4.  Velocities  of  ground  motion  at  elleven  different  sites  for 
the  event  Chancellor.  The  records  have  been  scaled  by 
epicentral  distance. 
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Figure  5.  Estimates  of  the  isotropic  part  of  the  moment 
tensor  for  the  events  Harzer  (top)  and 
Chancellor  (bottom). 
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DETERMINATION  OF  SOURCE  FUNCTION  FOR  A  SEISMIC  SWARM 


Pierre  MECHLER  -  Florence  RIVIERE 
RAD I OMAN A 


In  despite  of  the  title  of  this  presentation,  we  will 
not  present  a  study  of  source  functions,  but  a  way  to  obtain 
a  signal  averaged  over  the  seismic  stations  of  the  French 
Seismic  Network,  which  will  enable  us  to  study  the  source  in 
a  future  work. 


I .  Obtention  of  the  averaged  signal. 

A  classical  way  of  averaging  signals  over  a  seismic 
network  is  to  use  a  "Delay  and  Sum".  In  this  procedure,  a 
simple  stack  of  the  various  signals  are  done,  without  taking 
into  account  the  change  in  the  records  of  each  station. To 
overcome  this  problem,  we  decide  to  define  what  we  will  call 
the  best  averaged  signal  in  a  least  square  sens. 

The  technic  used  is  far  from  beeing  new.  We  suppose 
each  of  the  records  to  consist  of: 

-  a  common  signal  but  multiplied  by  a  station's  factor. 

-  a  noise  (mostly  wave  generated)  completely  decorrelated 

from  the  common  signal. 

The  signals  are  first  s’  'possed  to  be  delayed,  so  to  be 
in  phase,  but  we  will  cc^e  back  later  to  this  point.  We 
compute  the  correlation  matrix  of  all  signals.  It  can  be 
shown  that  the  station's  factors  is  one  of  the  eigen  vectors 
of  this  matrix,  the  one  corresponding  to  the  largest  of  all 
eigenvalues.  This  introduce  an  undetermination:  an  eigen 
vector  is  allways  obtained  up  to  a  constant,  but  this  is 
obvious  on  physical  basis,  it  is  possible  to  multiply  the 
station's  factors  by  an  arbitrary  number,  if  in  the  same 
time,  the  common  signal  is  divided  by  the  same  value. 

The  common  signal  is  then,  the  sum  of  the  records, 
averaged  by  the  just  computed  station's  factors.  The  wave 
generated  noise,  or  residual  signal  in  a  given  station,  is 
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obtained  by  a  difference  between  the  record  in  the  station 
and  the  common  signal. 

This  decomposition  of  the  signals  is  quite  similar  to 
the  covariance  analysis,  well-known  in  statistical  studies. 
The  main  difference  is  that  in  statistics,  the  decomposition 
of  the  signals  is  usually  done  on  a  larger  set  of  eigen 
vectors  than  just  the  first  one.  But,  in  our  case,  only  the 
first  as  a  meaning,  corresponding  to  a  signal  beneath  France 
and  representing  the  source  function  filtered  by  the 
propagation.  A  decomposition  on  a  larger  set  will  have  no 
physical  meaning,  the  coda  in  each  station,  are  only 
statisticaly  similar.  Of  course,  such  a  decomposition  has  to 
be  done  over  a  seimic  network  which  aperture  is  small  enough 
compared  to  the  epicentral  distance. 

This  procedure  is  relatively  sensitive  to  the 
determination  of  the  arrival  time  in  each  station.  We  were 
greatly  helped  in  this  by  the  origin  of  the  events  we 
studied.  It  consist  in  a  '’seismic  swarm”,  a  serie  of  very 
similar  events  from  the  same  area  in  Eastern  Kazakh,  not  far 
from  Shagan  River.  The  signals  recorded  in  one  station  for 
the  various  events  are  more  similar  than  the  signals 
recorded  in  the  stations  of  the  French  Seismic  Network  for 
the  same  event.  So,  in  a  first  step,  we  computed  the  common 
signal  in  each  single  stations  for  a  set  of  10  events.  It 
was  easier  to  measure  the  arrival  time  on  this  averaged 
signal,  then  to  obtain  a  precise  arrival  time  for  each 
record  by  a  correlation  technic.  We  also  checked,  by 
iteration,  that  it  was  not  possible  to  improve  the  time 
determination . 

We  vary  the  length  of  the  time  window  used  in  the 
computation.  The  samplig  rate  of  our  records  is  50Hz,  with 
an  anti-aliasing  filter  at  16Hz,  in  fact  all  of  our  data 
were  low-pass  filtered  at  8Hz.  The  various  time  windows  were 
between  1024  points  (about  20  seconds)  and  128  points  (2.5 
seconds)  without  real  change  in  the  station's  factors. 


II .  Magnitude  determination. 

The  first  interpretation  is  to  consider  the  common 
signal,  in  the  various  stations,  as  representing  the  signal 
arriving  beneath  France  and  filtered  by  an  average  crust, 
the  station's  factors  as  a  measure  of  the  sensitivity  of 
each  station,  and  the  residual  signal  as  athe  coda  (in  fact 
we  obtain  more  than  the  coda,  we  also  have  energy  in  the 
beginning  of  the  signal,  see  later). 
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We  will  start  by  magnitude  considerations.  We  studied 
two  different  ways  to  measure  the  magnitude.  In  the  first, 
for  each  station,  we  computed  for  all  events  a  common  signal 
and  ,in  that  case,  the  event's  factors.  The  magnitude  beeing 
the  logarithm  of  the  wave  amplitude,  plus  a  distance 
correction,  and  as  all  events  are  practicaly  at  the  same 
distance,  the  logarithm  of  the  event's  factors  should  be, 
within  a  constant  correction,  the  event's  magnitude. 

Using  the  magnitude  determination  of  the  French 
Network,  we  found  the  best  relation  between  magnitude  and 
log  of  event's  factor  to  be  of  the  form: 

mb  =  a*log(EF)  +  b 

but  with  a  "a"  factor  not  equal  to  1.  We  give  in  annex  the 
list  of  all  the  results,  the  mean  values  are: 

a  *  0.837  +/- . 086 

b  ■  5.952  +/- . 035 


We  also  used  only  nine  of  the  events  to  recompute  a  and 
b  ,  to  calculate  the  magnitude  of  the  last  one.  Volontarily, 
we  select  the  nine  smallest  events  to  obtain  the  magnitude 
of  the  largest.  The  result  were: 

m  =  6.17  +/— . 03 

which  is  exactly  the  value  found  by  the  French  Network. 

We  intent  to  apply  the  same  technic,  on  a  much  larger 
data  set,  in  order  to  obtain,  not  really  a  magnitude,  but  a 
classification  of  the  various  events,  according  to  their 
energy. 

This  first  way  to  determine  a  magnitude  is  rather 
tedious:  each  time  a  new  event  is  aded  to  the  data  base,  all 
computations  have  to  be  redone.  An  other  way,  is  to  use  the 
mean  station's  factors.  We  noticed  that  the  station'  factors 
are  very  similar  for  a  given  station  and  the  various  events. 
So,  we  first  computed  an  average  value  of  the  station's 
factor,  in  each  station  for  all  events,  then  we  computed  an 
approximation  of  the  common  signal,  using  this  new  set  of 
station's  factors  and  measured  it's  amplitude.  The  results 
are  given  in  annex. 
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III.  Source  function. 


The  common  signal,  from  all  stations  of  the  network  and 
a  given  event,  represents  an  approximation  of  the  incoming 
signal,  beneath  France  and  filtered  by  the  mean  french 
crust.  We  present  here  some  of  the  common  signal  we 
obtained.  In  this  preliminary  part  of  the  work,  we  have  not 
yet  compare  our  results  with  a  synthetic  sismogramm,  using  a 
source  function  and  a  propagation  effect. 

In  the  following  figures,  we  plotted  the  common  signals 
for  all  ten  events  used  in  this  study,  and  also,  for 
comparison,  the  actual  signals  in  some  stations  both  for 
the  best  event  and  an  average  one. 


IV.  Signal  residuals. 

The  difference  between  the  actual  record  and  the  common 
signal  represents  more  than  the  coda.  As  we  pointed  earlier, 
we  also  have  energy  in  the  first  part  of  the  signal. 

In  the  figures,  we  show  the  mean  signal  computed  using 
all  events  in  the  fives  stations  of  Massif  Central,  in  the 
middle  of  France.  And,  for  the  same  stations,  the  mean 
residual  signal. 

Ci  this  last  plot,  some  coherency  between  the  signals 
accross  the  subarray  is  clearly  seen. 

We  think  that  the  size  of  the  seismic  network  is  an 
important  parameter  for  the  coherency  of  the  waves  accross 
it.  If  the  aperture  is  too  small,  we  expect  to  obtain  a  high 
coherency  for  the  signals  themselves  (up  to  the  amplitude 
which  is  known  to  varry) ,  so  the  residual  signals  should 
represent  the  coda.  If  the  aperture  is  larger,  then  it  is 
possible  to  observe  propagations  through  the  array.  We  know, 
from  other  studies,  that  the  crust,  just  south  of  the  sub¬ 
array,  presents  larg  anomalies.  We  whish  to  correlate  both 
phenomens. 


V.  Preliminary  conclusions. 

This  presentation  is  only  a  presentation  of  a  current 
work.  Our  method  for  precise  magnitude  determination  has 
already  given  us  very  good  preliminary  results,  and  we  will 
improve  our  statistic  by  increasing  the  number  of  events. 
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The  common  signals,  obtained  for  a  given  quake,  over 
our  network  will  allow  a  study  of  the  source  function.  And 
we  really  hope  to  explain  the  main  trends  of  the  residual 
signals  by  correlation  with  the  geological  structure  around 
the  network. 
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Annex 


Magnitude 

a 

0.7527 

0.8390 

0.8184 

0.7185 

0.8635 

0.7973 

0.8387 

0.8313 

0.9131 

0.7698 

0.8024 

0.9291 

0.8525 

0.8399 

0.8257 

0.8985 

0.8849 

0.8236 

0.8575 

0.8890 

0.8103 

0.8578 

0.8799 

0.8114 

0.7131 

1.1354 

0.6380 


coef.  with  all  events 
b 

5.9658 

5.9690 

5.9445 

5.9974 

5.9124 

5.9924 

5.9021 

5.9924 

5.9124 

5.9687 

5.9529 

5.9509 

5.9552 

5.9574 

5.9629 

5.9497 

5.9441 

5.9285 

5.9212 

5.9227 

5.9494 

5.9210 

5.9135 

5.9365 

5.9693 

5.9337 

6.0732 


Table  I 


57 


Magnitude  computed  from  common  signal 


m=a*log (A) +b 


a — 0 . 886 

29/03/81 

13/09/81 

31/08/82 

26/12/82 

20/11/83 

06/10/83 

07/03/84 

10/02/85 

29/03/84 

14/07/84 


b=6 .406 

5.49 

6.20 

5.34 

5.60 

5.36 

6.03 

5.52 

5.69 

6.00 

6.18 


Table  II 


Magnitude  computed  from  station's  factors 


station 

eO 

el 

e2 

e3 

e4 

e5 

e6 

e7 

e8 

e9 

LBF 

5.56 

6.15 

5.43 

5.54 

5.24 

6.03 

5.40 

5.88 

5.97 

6.17 

LOR 

5.50 

6.14 

/ 

5.62 

5.35 

6.03 

5.42 

5.92 

5.97 

6.17 

AVF 

5.66 

6.13 

5.17 

5.57 

5.37 

5.87 

5.61 

5.85 

6.02 

6.18 

LSF 

/ 

/ 

/ 

/ 

5.33 

/ 

5.48 

5.88 

6.01 

6.19 

SSF 

5.60 

6.12 

5.09 

5.63 

5.43 

5.91 

5.65 

5.90 

5.97 

6.12 

MFF 

/ 

6.15 

5.44 

5.50 

5.26 

6.02 

/ 

5.91 

5.98 

6.17 

GRR 

/ 

/ 

5.17 

5.63 

5.39 

5.92 

5.62 

/ 

5.97 

6.15 

FLN 

/ 

6.20 

5.31 

5.58 

5.33 

5.94 

5.47 

/ 

6.02 

6.17 

LDF 

/ 

/ 

5.24 

5.61 

5.37 

5.95 

5.51 

/ 

6.01 

6.17 

SMF 

5.63 

6.19 

5.39 

5.52 

5.28 

/ 

5.42 

5.87 

5.98 

6.17 

LRG 

5.67 

6.16 

5.20 

5.57 

5.33 

5.85 

5.59 

5.85 

6.02 

6.18 

TCF 

/ 

6.17 

5.25 

5.64 

5.29 

5.98 

5.55 

5.93 

5.98 

6.13 

RJF 

5.63 

6.18 

5.26 

5.54 

5.34 

5.92 

5.48 

5.87 

6.02 

6.18 

FRF 

5.63 

6.17 

5.30 

5.53 

5.31 

5.93 

5.47 

5.91 

5.99 

6.19 

LPO 

5.61 

6.17 

5.30 

5.53 

5.34 

5.95 

5.45 

5.87 

6.00 

6.19 

LFF 

5.58 

6.18 

5.26 

5.56 

5.34 

5.96 

5.48 

5.88 

5.99 

6.18 

EPF 

5.59 

6.16 

5.24 

5.57 

5.36 

6.05 

5.49 

5.82 

5.96 

6.17 

BSF 

5.60 

6.11 

5.14 

5.59 

5.38 

5.96 

5.58 

5.93 

6.00 

6.13 

HAD 

5.62 

6.16 

5.20 

5.58 

5.36 

/ 

5.50 

/ 

5.98 

6.16 

CDF 

5.55 

6.17 

5.15 

5.65 

5.37 

/ 

5.60 

5.90 

5.97 

6.10 

LMR 

5.64 

6.13 

5.22 

5.53 

5.31 

5.94 

5.57 

/ 

6.04 

6.15 

LPF 

/ 

/ 

5.23 

5.59 

5.37 

5.92 

5.56 

5.87 

6.02 

6.19 

CAF 

/ 

/ 

5.25 

5.62 

5.34 

5.99 

5.50 

/ 

5.99 

6.16 

MZF 

/ 

/ 

5.30 

5.53 

5.34 

5.93 

5.53 

5.89 

6.03 

6.19 

CVF 

/ 

/ 

/ 

5.57 

/ 

5.96 

5.54 

5.88 

6.01 

6.18 

BGF 

/ 

/ 

/ 

/ 

5.27 

/ 

5.55 

5.91 

5.98 

6.17 

HYF 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

5.90 

6.00 

6.18 

mean 

5.60 

6.16 

5.24 

5.58 

5.34 

5.95 

5.52 

5.89 

5.99 

6.17 

sig 

0.04 

0.02 

0.09 

0.04 

0.04 

0.05 

0.06 

0.02 

0.02 

0.02 

LDG 

5.5 

6.18 

5.3 

5.6 

5.3 

5.96 

5.5 

5.89 

6.02 

6.17 

Table  III 
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Figures  Caption 


Fig  1  :  common  signals  for  the  ten  events  of  this  study 

Fig  2  :  Actual  signals  for  event  1  :  13/09/81 

Fig  3  :  Actual  signals  for  event  4  :  20/11/83 

Fig  4  :  Averaged  records  in  the  5  stations  of  Massif  Central 

Fig  5  :  Averaged  residual  signals  in  the  5  stations  of 

Massif  Central 


60 


event  9 


event  o 


event  7 


I 

i 


event  6 


event  d 


!  I 

M 

I1, 

event  4  ^  [  \ 

"  — :i  I  ‘  v*"’ 

y . 

event  3  i i 


v , 


event  2  p, 

- -~-,n  /w 


event  1 


•  :  i  r 
\  »  .1 


event  0  (\ 

_ _  >i 


Fig  1 


61 


62 


LMR 

BSF 

EPF 

LFF 

LPO 

FRF 

RJF 

MFF 

SFF 

LOR 


I 

rig  3 


63 


Fig  4 


64 


SMF 


i 


I1' 

.‘i  i 

i  <  -1  i 

!  Ml  I 


SSF 


A  VF 


LOR 


Fig  5 


65 


THE  DAMAGE  MECHANICS  OF  BRITTLE  SOLIDS  IN  COMPRESSION 


C.  G.  SAMMIS1  and  M.  F.  ASHBY2 


department  of  Geological  Sciences 
University  of  Southern  California, 

Los  Angeles,  CA  90089-0741,  USA  and 
department  of  Engineering,  University  of  Cambridge, 
Trumpington  Street,  Cambridge  CB2  IPZ,  U.K. 


66 


INTRODUCTION 


When  a  brittle,  Inhomogeneous  solid  is  loaded  in  compression,  microcracks 
nucleate  and  propagate  from  the  inhomogeneities.  The  difference  between 
compressive  and  tensile  fracture,  which  is  much  better  understood,  is  that  a 
brittle  tensil  crack  is  unstable  (one  crack,  once  started,  accelerates  across 
the  sample  to  give  failure);  by  contrast,  brittle  cracking  in  compression  is 
at  first  stable,  so  that  many  small  cracks  extend,  each  growing  longer  as  the 
stress  is  raised,  until  they  interact  in  some  co-operative  way,  to  give  final 
failure . 

Systematic  measurements  of  strength  really  began  about  the  middle  of  the 
last  century  (see,  for  example,  Jaeger  and  Cook,  1976)  but  without  much 
attempt  to  understand  what  determined  it,  or  why  the  compressive  strength  is 
ten  or  more  times  greater  than  that  in  tension.  Progress  in  understanding 
these  problems  is  more  recent.  A  series  if  recent  papers  and  reviews  (Griggs 
and  Handin,  1960;  Paterson,  1978;  Hallbauer  et  al.,  1973;  Tapponnier  and 
Brace,  1976;  Wawersik  and  Fairhurst,  1970;  Wawersik  and  Brace,  1971; 
Nemat-Nasser  and  Hori,  1982;  Kobayashi,  1971;  Newman,  1978;  Ashby  and  Hallam, 
1986;  Sammis  and  Ashby,  1986)  have  established  that,  in  compression,  numerous 
small  cracks  nucleate  and  grow  stably  (meaning  that  each  increment  of  crack 
growth  requires  an  increment  in  load  to  drive  it).  An  isolated  crack  in  a 
very  large  body  continues  to  grow  in  this  stable  way  until  the  length  becomes 
comparable  with  that  of  the  body  itself,  but  the  multitude  of  cracks  which 
nucleate  in  most  natural  rocks  behave  differently:  when  the  crack  length  is 
comparable  with  their  spacing  they  interact  and  link,  precipitating  failure. 
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The  problem  is  complicated  by  time-dependent  effects  (Anderson  and  Grew, 
1976;  Martin,  1971;  Waza  et  al.,  1980,  Sano  et  al.,  1981,  Costin  and  Holcomb, 
1981).  They  have  at  least  two  origins.  First,  cracking  in  compression  is 
associated  with  dilation;  if  the  body  is  saturated  with  a  fluid,  then  its  flow 
into  the  dilating  region  can  introduce  a  time-dependent  aspect  to  fracture. 
Second,  in  some  brittle  solids,  notably  siliciates,  crack  growth  is  slow, 
limited  by  chemical  reaction  (often  with  water)  at  the  tips  of  the  cracks,  as 
well  as  by  the  rate  of  loading.  In  both  cases,  a  static  load  which  does  not 
immediately  cause  failure  may  still  do  so  if  left  in  place  for  a  sufficient 
length  of  time. 

The  understanding  of  these  phenomena  is  still  incomplete.  But  the 
physical  basis  of  compressive  brittle  fracture  is  much  clearer  than  it  was  a 
decade  ago.  It  seems  an  appropriate  time  to  tty  to  abstract,  from  the 
obser/ations  and  modelling  which  are  currently  available,  a  simplified 
descri  >tion  of  compression-cracking,  retaining  as  far  as  possible  the 
important  physical  understanding.  The  goal  is  to  develop  a  damage  mechanics 
of  brittle  solids,  from  which  the  stress-strain  response  and  an  operational 
defin'cion  of  failure  can  be  derived  for  any  given  stress  state  and  initial 
set  o*  material  properties.  Two  attempts  to  achieve  this  can  be  found  in  the 
open  ’iterature;  that  of  Costin  (1983),  and  that  of  Sammis  and  Ashby  (1986). 
Central  to  the  problem  is  the  relationship  between  stress  and  crack  extension. 
Costin  (1983)  postulates  an  assumed  relationship  of  reasonable  form,  and 
develops  from  it  expressions  for  the  filure  surface  which  (with  some 
adjustable  parameters)  give  a  good  description  of  a  body  of  experimental  data. 
Sammis  and  Ashby  (1986)  and  Ashby  and  Hallam  (1986)  developed  a  detailed 
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physical  model  for  crack  extension,  which  they  use  to  plot  stress-strain 
curves  for  brittle  solids  from  which  failure  surfaces  can  be  constructed;  but 
the  complexity  of  the  model-based  equations  makes  the  process  cumbersome.  In 
the  present  paper  we  attempt  to  develop  a  simpler,  model-based,  mechanics  of 
brittle  compressive  fracture,  drawing  heavily  on  these  two  previous  pieces  of 
work. 

A  model-based  damage  mechanics  has  the  advantage  that  several  important 
phenomena  can  be  incorporated  in  a  natural  way.  The  failure  criterion  can  be 
defined  in  terms  of  the  peak-stress  in  the  stress-strain  curve  rather  than 
requiring  the  postulation  of  a  critical  damage  as  in  Costin  (1983).  The 
initial  porosity  is  included  in  the  modeling,  as  is  the  size  spectrum  of  the 
intial  defects.  Finally,  the  effect  of  pore  fluids  on  the  evolution  of 
microcrack  damage  can  also  be  incorporated. 

We  have  not  yet  formulated  a  complete  damage  mechanics,  nor  have  we 
tested  is  against  the  wealth  of  compressive  failure  data  in  the  literature. 
This  paper  is  thus  a  progress  report;  but  it  contains  enough  detail  to  show 
the  strategy  and  framework  of  the  approach  we  are  using. 

THE  CRACK  INITIATION  SURFACE 

Almost  all  brittle  solids  contain  inhomogeneties :  holes,  little  cracks, 
particles  which  are  poorly  bonded  or  which  have  different  moduli  from  the 
matrix.  Any  of  these  can  act  as  nuclei  for  new  cracks  when  the  solid  is 
loaded. 

The  range  of  possible  nuclei  is  wide,  but  the  spectrum  of  their 
characteristics  is  probably  bracketted,  at  least  approximately,  by  two 
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extremes:  the  spherical  hole,  and  the  sharp  inclined  crack  (Fig.  1).  Both 

have  been  studied  experimentally  and  modelled,  the  first  by  Sammis  and  Ashby 
(1986),  the  second  by  Nemat-Nasser  and  Horii  (1982)  and  Ashby  and  Hallam 
(1986),  and  explicit  equations  for  the  initiation  of  cracking  are  developed 
for  both.  Details  of  the  results  are  given  in  Appendix  1;  here  we  note  simply 
that  both  lead  to  results  which  are  almost  identical,  and  can  be  written,  for 
axisymmetric  loading,  as  a  linear  relationship  between  the  axial  stress  o^  and 
the  radial  stress  02  =  03-  For  initiation: 

al  =  C1C3  +  C'l  (!) 

When  cracks  nucleate  from  spherical  holes,  the  constant  (Appendix  1) 
is  predicted  to  be  3.4.  When  they  nucleate  from  cracks,  Cx  depends  on  the 
coefficient  of  friction  between  the  sliding  crack  faces;  taking  p=0.6  as 
typical  for  rocks  gives  Cx  =  3.1. 

The  constant  C2  depends  principally  on  the  fracture  toughness  Kxc  of  the 
brittle  solid,  and  on  the  size  scale  of  the  inhomogeneity  (the  hole  or  crack). 
When  the  nuclei  are  holes,  C2  =  1.6  KIC/ /ir£  where  l  is  the  length  of  the 
incipient  cracks  associated  with  the  hole  (Fig.  1).  When  nuclei  are  inclined 
cracks,  C2  =  3.1  Kxc/*^.  where  2a  is  the  inclined  crack  length  (Fig.  1). 
Despite  the  rather  different  geometry,  it  is  obvious  that  not  only  is  the  form 
of  the  nucleation  equation  identical,  but  even  the  constants  it  contains  are 
almost  the  same  for  holes  and  cracks,  encouraging  the  view  that  Equation  (1) 
has  generality  as  an  equation  for  damage  initiation. 

Crack  initiation  can  be  detected  in  several  ways:  by  the  start  of 
acoustic  emission,  by  the  first  non-linearity  of  the  stress-strain  curve,  or 
by  a  sudden  increase  in  internal  friction.  Figure  2  shows  data  for  crack 
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initiation  obtained  by  the  first  two  techniques  for  Westerly  granite.  Brace, 
et  li* »  1966)  plotted  on  axes  of  03  and  03  to  allow  comparison  with  Equation 
(1).  The  linear  relationship  is  clearly  a  good  description  of  the  data.  The 
slopes  C3 ,  in  both  sets  of  data  are  close  to  4.1  (the  models  suggest 
3.2  -  3.4),  and  the  constants  C2,  in  both  cases,  give  l  »  a  *  0.6  mm  (taking 
Kjc  =  1  MPa  mV2),  which  is  reasonable. 

THE  EVOLUTION  OF  DAMAGE 

As  the  body  containing  pores,  angled  microcracks,  etc.,  is  loaded,  the 
stress  intensity  at  the  tips  of  the  wing-cracks  grows.  We  will  assume  that 
the  wing-cracks  advance  so  that,  at  their  tips,  Kj  =  Kjc>  though  in  reality 
the  inhomogeneity  of  the  material  may  cause  them  to  extend  in  little  jumps. 

The  problem,  then,  is  to  estimate  Kx  at  the  tip  of  the  wing  cracks. 

Crack  Growth  and  Interation 

The  development  of  a  stress  intensity  at  the  wing-crack  tip  caused  by  a 
stress  field  can  be  likened  to  the  induced  polarisation  of  an  insulator  by  an 
electric  field.  First,  there  is  the  polarisation  of  an  isolated  atom;  but  in 
addition,  the  polarisation  of  neighboring  atoms  adds  to  the  field  and  causes 
further  polarisation  of  the  first  atom.  In  a  rather  similar  way,  the  remote 
stress  field  (03,  02,  03)  induces  a  contribution  to  Kx  (which  we  shall  call 

OO 

Kx)>  But  the  cracks  perturb  the  stress  field  inside  the  body,  so  that  a  given 

i 

crack  experiences  an  additional  contribution  to  Kx  (which  we  call  Kx)  from  its 
interaction  with  the  stress  fields  induced  at  its  neighbors.  The  problem  is 
like  that  of  induced  polarisation  because  the  internal  stresses  only  appear  in 
response  to  the  external  loads. 
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Consider  first  the  stress  intensity  Kj  which  is  directly  induced  at  an 
isolated  crack  in  an  infinite  medium  by  the  applied  stress.  Its  value  depends 
on  the  geometry  of  the  defect  from  which  the  wing  cracks  grow.  When  this  is 
an  angled  crack,  the  result  for  a  crack-length  L  »  1  is  (Ashby  and  Hallam, 


1986,  eqn.  (6)): 


„  C3  01  /  -rra 

Ki  =  -  - -  +  03/ ira  L  1'2 

1  L  1/2  +  3 


Coal'S  va  XL 

— r  - 1 


with  C3  =  0.23  and  L=£/a.  Figure  20  of  Ashby  and  Hallam  (1986)  and  Figure  20 
of  Sammis  and  Ashby  (1986)  both  show  that  the  peak  stress  occurs  at  3<L<5,  so 
the  approximation  L»1  is  tolerable.  The  first  of  the  two  terms  describes  the 
stress  intensity  caused  by  the  axial  stress  a\\  it  decreases  as  the  crack 
extends,  giving  stable  crack  growth.  The  second  term  describes  the  closing 
effect  of  the  lateral  stress  03;  it  acts  on  the  whole  crack  length  L  and  thus 
grows  in  importance  as  L  increases.  The  analagous  result  for  a  crack  growing 
from  a  hole  (Sammis  and  Ashby,  1986,  eqn.  (4)  and  (5))  is  a  little  more 
complicated.  In  the  limit  of  L  »1,  we  have: 


C4  ai  /rra 


+  03/7?^  L  1^2 


C4  0i/ira  1  -  XL 
L  1/2  L2  C4  ^ 

with  C4  ■  1.1.  The  interpretation  of  the  two  terms  is  the  same  as  before. 
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We  now  consider  the  interaction  between  cracks,  giving  the  additional 
i 

contribution  Kj.  In  our  earlier  work  we  modelled  interaction  by  considering 
the  bending  displacements  between  neighboring  cracks.  This  idealization  was 
particularly  successful  in  modelling  the  interaction  of  a  crack  with  a  freee 
surface,  but  it  is  realistic  for  the  interaciton  between  cracks  only  when  they 
overlap  extensively.  To  avoid  this  assumption  we  develop  a  new  model,  as 
follows . 

Consider  a  starter  crack  or  void  with  vertical  wing  cracks  of  length  £. 

When  the  axial  stress  is  applied,  sliding  (at  an  angled  crack)  or  bulging 

(at  a  void)  wedges  the  wing-cracks  open,  that  is,  it  induces  stress  intensity 
00 

Kj  at  the  wing  crack  tips.  This  crack  is  now  a  source  of  an  elastic  field  in 
the  body.  We  measure  the  induced  strength  of  this  crack  by  the  mid-point 
crack  opening  displacement,  6.  The  crack  is  like  a  dislocation  loop  of  radius 
£  and  Burger's  vector  b  =  5,  (though  note  that  6  is  induced  by  the  external 
loads  and  thus  increases  with  stress , whereas  the  strength  of  a  dislocation 
loop  is  fixed).  The  stresses  at  a  distance  r  »  £  from  a  dislocation  loop  or 
dip  2e  scale  as: 

GbZ 

a  =  C  -  (4) 

r2 

Our  model  is  that  the  locally  induced  stresses  of  a  crack  of  radius  £  and 
strength  5,  at  a  distance  r  »  £  from  the  crack  center  scale  as: 

ES£ 

ac  =  C5  -  (5) 

(r-£) 2 

where  we  might  expect  C5  *  0.5, 

The  stress  obviously  varies  with  angle  around  the  crack;  in  some  places 
it  will  have  components  which  will  tend  to  open  other  cracks,  in  some  places 
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it  will  tend  to  close  them.  We  focus  on  the  region  in  which  the  crack  field 

has  an  opening  component.  A  second  crack  lying  in  this  region  now  experiences 
00  i 

a  Kj  from  the  applied  stress  and  a  Kj  from  the  induced  field  of  the  first 
crack.  Its  magnitude  is: 


i  _ 

Ki  =  Oc/tt£  (6) 

It  remains  to  evaluate  ac,  insert  the  result  into  this  equation,  and  add  the 
two  contributions  to  Kj . 

The  crack  opening  displacement  6  is  known.  Both  when  the  wing  cracks 

grow  from  voids  and  when  they  grow  from  angled  cracks  (or  from  anything  else) , 
TOT 

5  is  related  to  Kj  by: 


Kitot/£ 

5  =>  eg  -  (7) 

E 

when  cf,  =  4  (Tada,  Paris  and  Irwin,  1973).  (This  expression  is  completely 
general,  and  does  not  depend  on  details  of  crack  origins  or  loading).  Thus, 
combining  Equations  (4),  (5)  and  (6) 

£2 

Ki101  =  Ki"  +  Kx1  =  Ki'Va  -  C7  -  2  )  (8) 

(r-£) 

The  second  term  in  the  brackets  describes  the  interaction,  and  Cj  combines 
the  constants  C5 ,  C^,  /x  etc.  It's  value  is  roughly  C7  “3.5. 

The  total  stress  intensity  at  a  "typical"  crack  can  now  be  written  in 
dimensionless  variables  L  =  £/a  and  R  =  r/a  as  follows.  For  wing  crack  growth 
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from  an  angled  crack: 


Inverting  these  equations  gives,  for  angle  cracks: 


iv  iv 

L  L  l/2 

L2 

°1  =  °3 

C3  C3 

(l-c7 

(R-L)2 

and,  for  holes: 

L3  L  5/2 

L2 

°1  =  °3 

c4  c4 

(1-C7 

(R-L) 2 

(9) 


(10) 


(ID 


(12) 


where  a  =  o/xa/Kjc  is  the  dimensionless  stress.  For  holes,  ira2/r2  is  the  area 
of  holes  per  unit  area  of  sample,  or  the  porosity.  Denoting  the  porosity  by 
fa,  we  thus  have  R  *  (w/fg)1/2.  When  fractures  are  initiated  at  angle  cracks, 
fa  is  interpreted  as  that  fraction  of  the  total  area  which  contains  starter 
cracks,  where  an  area  7ra2  is  associated  with  each  starter  crack  of  length  2a. 
Again  R  =  (ir/fa)  ^ 2 • 

It  is  convenient  to  define  damage  0<  D  <  l  as 

D=£2Na=L2fa/ir 

Where  Na  is  the  number  of  starter  flaws  per  unit  area. 

Equations  (11)  and  (12)  are  plotted  In  Figure  3  for  two  different 
confining  stresses  at  fixed  for  initial  porosity  (or  crack  density).  The 
stress  passes  through  a  maximum  and  then  falls  as  the  damage  increases. 
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The  peak  stress  corresponds  to  one  definition  of  failure.  The  failure 
envelope  for  Westerly  granite  in  Figure  2  was  generated  by  plottin6  peak  03, 
as  a  function  of  03. 

Damage  also  lowers  the  modulus  of  the  cracked  body.  If  the  modulus  of 
the  solid  is  E0,  then  introducing  the  holes  or  initial,  inclined  cracks,  lower 
the  modulus  to  E  =  E0/(l+2fa),  where  fa  =  Naira2  is  the  area  fraction  of  holes 
or  cracks.  The  cracking  which  occurs  on  loading  further  lowers  the  modulus. 
Sammis  and  Ashby  (1985)  show  that  the  modulus  may  be  approximated 

Eo 


E  = 


l+2fa+2D  (1-3 X) 2 


(13) 


The  strain  is  simply  e  =  a/ E,  or 


ai 

e  =  —  (l+2fa+2D  (1-3 X) 2)  (14) 

E0 

These  results  can  be  used  to  calculate  the  stress-strain  curve,  and  the 
failure  strain. 


PLASTIC  YIELDING 

As  the  confining  pressure  is  increased,  brittle  fracture  is  made 
increasingly  difficult.  A  critical  pressure  may  be  reached  at  w >ich  true 
plasticity  replaces  crack  extension.  This  transition  can  be  illustrated  by 
plotting  a  yield  (or  a  creep)  surface,  defined  by: 

1 

-  ((03  -  c^)2  +  (CT2  ~  CT3)2  +  ( a3  ~  al)2)  =  ay2  (15) 

2 

This  is  shown  as  broken  lines  on  Figure  2.  The  yield  strength,  Oy,  can  be 
derived  from  hardness  data,  since  Oy  *  h/3. 
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DISCUSSION 


As  pointed  out  in  the  introduction,  our  damage  mechanics  formulation  is 
not  yet  complete  or  tested.  We  conclude  this  progress  report  by  discussing 
some  areas  which  require  additional  work. 

Crack  Growth  from  a  Distribution  of  Nuclei  of  Different  Sizes 

In  most  solids,  the  starter  flaws  are  not  all  the  same  size.  This  is 
evident  from  the  observed  change  in  failure  mode  associated  with  increasing 
confirming  pressure.  For  unconfined  samples,  a  few  of  the  largest  flaws  are 
activated  and  grow  until  they  intersect  the  sample  ends.  The  sample  fails  by 
"vertical  slabbing".  The  application  of  even  a  small  confining  stress 
supresses  the  growth  from  these  larger  flaws,  and  allows  the  nucleation  of 
cracks  from  smaller  or  less  favorably  oriented  starter  flaws.  Before  any 
single  crack  can  grow  sufficiently  long  to  intersect  the  sample  ends,  the 
density  of  activated  flaws  has  increased  enough  to  produce  a  shear  instability 
and  the  sample  fails  by  faulting.  At  very  large  confining  pressures  the 
extreme  supression  of  crack  growth  allows  the  activation  of  so  many  flaws  that 
many  shear  instabilities  are  produced,  resulting  In  pseudo-ductile 
deformation. 

Another  indication  of  the  importance  of  including  a  range  of  flaw'-sizes 
in  our  model  is  the  negative  curvature  of  the  failure  surface  in  Figure  2. 

Note  t'.at  our  model  surface  is  approximately  linear  when  all  flaws  are  the 
same  size.  Activation  of  additional,  smaller  flaw  with  increasing  stress  will 
produce  the  observed  negative  curvature. 
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In  order  to  introduce  this  behavior  into  our  model,  we  will  assume  that 
the  size  distribution  of  starter  flaws  may  be  approximated  by  the  Weibull 
distribution. 

tot  _  n 

fa(a)  =  fa  exP  “  (a/a)  (16) 

tot  _ 

where  fa  is  the  total  area  fraction  of  starter  defects  of  all  sizes  and  a 

is  a  measure  of  the  average  size  of  the  initial  flaws.  The  parameter  n  sets 

the  width  of  the  distribution:  small  n  corresponds  to  a  broad  distribution 

while  for  large  n  all  flaws  have  a  *  a. 

The  number  of  flaws  which  have  been  activated  depends  on  the  stress  state 
( ai  and  03).  Writing  the  initiation  condition  as 


where  o 


activiated . 


Increasing  at  constant  03  causes  smaller  flaws  to  be 
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Note  that  the  definition  of  damage  given  above  for  the  case  where  all 
flaws  are  the  same  size  is  still  valid 

ia*(  cri ,  03)  ra*(a1,  03) 

ir£2dNa  =  L2  Jjra2dNa  =  L2fa  (01,  03)  (20) 

since  L  is  a  constant  for  given  ( ,  03).  The  damage  mechanics  will  be  as 
developed  above  except  that  fa  is  now  a  function  of  the  stress  state. 

Crack  Initiation  And  Growth  In  A  Saturated  Porous  Medium 

When  pores  are  saturated,  the  stress  state,  and  thus  the  stress  intensity 
factor  is  different  than  that  given  in  Equation  (3).  Fluids  in  the  initial 
inclined  cracks  might  be  expected  to  lower  the  effective  coefficient  of 
friction.  We  are  now  developing  these  models  for  comparison  with  "effective 
stress  laws"  in  current  use. 

Crack  Interaction 

Central  to  the  entire  damage  mechanics  approach  to  brittle  failure  in 
compression  is  the  crack  interaction.  Sammis  and  Ashby  (1986)  and  Ashby  and 
Hallam  (1986)  point  out  that  physically,  this  interaction  takes  the  form  of 
mutual  growth  enhancement  of  parallel  neighboring  cracks.  How  well  this 
phenomena  is  represented  by  Equatios  (5)-(7)  remains  to  be  tested.  We  are 
currently  testing  this  hypothesis  in  two  ways:  (1)  by  generating  crack 
initiation  and  failure  envelopes  (like  Figure  2)  for  the  wide  range  of 
materials  for  which  triaxial  data  currently  exist,  and  (2),  by  performing 
model  experiments  on  two-and  three-dimensional  arrays  of  holes  and  cracks. 

The  fact  that  Costin  (1983)  got  promising  results  postulating  a  critical 
damage  failure  criterion  and  ignoring  interactions  suggests  that  the  final 
results  may  not  be  too  sensitive  to  the  exact  functional  form  assumed  for  the 


Interaction. 
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APPENDIX  1:  INITIATION  OF  CRACKING 


(a)  Crack  nucleation  at  a  hole 

Sammis  and  Ashby  (1986)  give  the  criterion  for  the  initaiton  of  cracking 
from  a  spherical  hole  as: 


oi/ira  -(l+L)4*1 


KIC 


✓L  {0.62  (1  -  1 . 8 X)  -  A(l+L)4-!} 


(Al) 


Where  is  the  axial  stress,  X=a3/oi  is  the  ratio  of  radial  (03)  to  axial 
(cri)  stress,  a  is  the  hole  radius,  L=l/a  is  the  normalized  crack  length  and  l 
is  the  absolute  value  of  the  crack  length.  Kjq  is  the  fracture  toughness  of 
the  material.  For  initation  L  is  small  so  that  1+L  ■  1. 

Replacing  X  by  03/ o\  and  rearranging  gives: 

1.6  Kjc 

ai  =  3.4  03  +  ( A2 ) 

/L  /na 

Note  that  the  second  term  on  the  right  is  a  constant  for  a  given 
most-dangerous  flaw  characterized  by  a  hole  of  radius  a  with  a  starter  crack 
of  length  L  =  l/a.  The  equation  predicts  a  linear  relationship  between  oi  and 
03  with  a  slope  of  3.4  (or  thereabouts)  independent  of  the  flaw  distribution. 
The  intercept  allows  /£  to  be  calculated. 


(b)  Crack  Nucleation  At  An  Inclined  Crack 

Nemat-Nasser  and  Hori  (1982)  and  Ashby  and  Hallam  (1986)  analyse  the 
nucleation  of  wing  cracks  from  an  initial  inclined  crack.  The  condition  can 
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be  written  as: 


o^/ira  -V3 

—  —  ■  -  ■  s  —  —  ■'  I  — —  .  I  I  .1.  —  .  ■ 

KIC  {(l-X)(l+V2)l/2  -  U+A)y} 


(A3) 


where  2a  is  the  length  of  the  initial  inclined  crack,  y  is  the  coefficient  of 
friction,  and  the  other  symbols  are  defined  above.  Replacing  X  by  03 /  03  and 


rearranging  gives: 


( l+y2)  l/  2  +  y  KIC 

oi  »  { - - - )  03  +  — r — 

(1+y2)  1//2  -  y  /rra 


(AA) 


As  before,  the  second  term  on  the  right  is  constant  for  a  given 
most-dangerous  flaw,  characterized  by  a  length  2a.  The  coefficient  of  03  has 
a  value  between  3  and  5,  depending  on  the  value  of  y  (between  0.5  and  1). 


smc357 
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A 


FIGURE  1  Starter-flaws  for  which  crack  nucleation  and  growth  nas  Deen 
analysed.  The  sharp  angle  crack  and  spherical  pore  are  considered 
end -members  of  the  spectrum  of  flaw  geometries  in  rock. 
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FIGURE  3  Axial  stress  as  a  function  of  damage  for  two  values  of  the 
confining  stress.  The  peak  value  of ^  ts  taKen  as  the  failure  stress.  The 
broken  curves  were  calculated  using  Equation  (11)  for  angle  cracks  while 
the  solid  curves  were  calculated  using  Equation  ( 1 2)  for  holes. 
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In  the  last  annual  report  we  have  described  a  method  (  Campillo  and  ai  3984) 
to  extract  from  the  Lg  waves  train  the  Q  factor  associated  to  the  propagation 
region  and  the  source  spectrum. 

We  have  applied  these  results  for  some  quakes  occuring  in  Remiremont 
(  Vosges)  in  December  84/  January  85  and  recorded  on  the  seismic  sta*  >ns 
of  the  LDG  network.  The  aim  of  this  study  was  to  evaluate  the  source 
amplitude  spectrum  for  the  main  shock  and  some  precursors  and  after 
shocks,  and  check  if  these  results  extracted  from  Lg  waves  in  the  far 
field  were  coherent  with  those  obtained  at  short  distance  from  the  epicenter. 

Seismic  data 

Among  the  28  SP  seismic  stations  of  the  LDG  network  only  8  having  recorded 
this  crisis  with  a  reasonable  signal  upon  noise  ratio  have  been  selected. 

In  the  Figure  8  are  represented  the  data  recorded  for  the  main  shock 
(  12/29/84  ;  Ml  =  4.  3  )  and  the  stations  versus  epicenter  repartition. 

All  quakes  are  superficial,  between  5  and  10  km  depth.  Seismic  data  are 
numerically  recorded  on  magnetic  tape  with  a  sampling  rate  of  the 
50  samples  /  sec  which  allows  to  process  signals  spectra  up  to  15  Hz. 
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Results 


After  shock  01.02.85  I^L  =  4.  0 

This  after  shock  is  recorded  at  HAU  station  (  A  =  21  km  ) 

without  saturation. 

After  deconvolution  of  the  seismograph  response,  the  amplitude  spectrum 
of  the  Sg  signal  is  computed. 

It  represents  an  amplitude  spectrum  of  the  source. 

A  source  displacement  model  of  the  form  : 


n  (f) 


a  ( o ) 


V1  +  <■ 


i. 

fc 


(  Boatwright  1978  ) 


is  computed  by  least  squares  to  fit  this  amplitude  spectrum  (  Fig.  9  a.  ) 
We  find  corner  frequency  fc  =3.6  Hz 

slope  ^  1 

The  seismic  moment  Mo  is  also  available  by  the  following  relation  ship 
(  Street  and  al  197  5  ;  Hermann  and  Kijko  1983) 

s'tt'  p  a  3  r0  n  (o) 

Mo  =  - 5 - 

R  84> 

rQ  =  1  km 

jQ  (o)  ;  displacement  source  spectrum  at  very  low  frequency 
jb  :  shear  wave  velocity  =  3.5  km/s 
^  '  density  2.2  g/cm^ 

:  coefficient  taking  into  account  both  radiation  pattern  ^  ] 

20 

We  find  Mo  =  2.8x10  dyne  .  cm 

We  also  compute  from  Lgj  ,  Lg  2  ,  Lg  3  waves  trains  recorded 
in  the  8  stations  situated  at  more  than  200  km  from  the  epicentral  area, 
a  mean  amplitude  source  spectrum.  For  this  computation  we  correct 
the  Lgi  t  Lg  2»  Lgj  ,  spectra  from 

-  elastic  and  anelastic  attenuation 

-  station  response 

-  seismograph  response 
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Figures  9b,  9  c  ,  9  d  represent  an  evaluation  of  the  o”  limits 

for  the  8  source  spectra  and  the  mean  amplitude  source  spectrum. 

(  continue  line  ). 

We  also  figure  the  source  spectrum  computed  from  HAU  Sg  waves 
(  doted  line  )  as  a  reference. 

As  a  result  ,  source  spectra  obtained  with  Lg  waves  in  the  far  field  are 
similar  to  the  one  obtained  in  the  near  field  (  HAU  )  namely.  Lg]  waves 
lead  to  a  source  spectrum  which  JTI  (o),  fc  and  slope  are  very  comparable 
to  those  of  HAU. 

Seismic  moments  computed  on  the  source  spectra  from  Lg  waves  and  the 
one  extracted  from  Sg  waves  at  HAU  are  of  the  same  order. 


-  Main  shock  12.29.84  Ml  s  4.  8 

Source  amplitude  spectra  on  Lg  ]  t  Lg2  ,  Lg  3  waves  are  computed 
by  the  same  method  described  in  the  first  case  (  After  shock  01  02  85 
ML  =  4.0  ). 

For  this  main  shock  case  no  records  are  available  at  HAU. 

Figure  10  represents  the  source  spectra  we  obtain  : 

-  -  C~  amplitude  spectra  corresponding  to  all  the  8  stations 

-  the  source  spectrum  which  fits  the  mean  amplitude  spectrum  values. 

-  corner  frequency  fc  ,  slope  of  the  spectrum  and  seismic  moments 
which  are  very  similar  . 

-  fault  dimensions  and  stress  drops  computed  for  each  of  the  three  cases 
with,  for  the  circular  fault  radius,  two  different  formula  proposed  by  Brune 
and  Madariaga. 


L 


0.  37 


l 

fc 


L 


O.  21 


(  Brune  ) 


(  Madariaga  ) 


and 


Mo 
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-  Small  foreshock  12.24.84  Ml  =  3.0 

This  last  example  of  source  evaluation  is  given  to  figure  out  the 

limits  of  our  far  field  method.  This  small  foreshock  is  recorded  by  8  stations  j 

among  those  which  recorded  the  previous  cases. 

Figure  11  presents  the  source  amplitude  spectra  obtained  at  HAU  on 
Sg  waves  and  on  Lg]  ,  Lg2  .  Lg3  waves  recorded  on  the  8  stations 
at  200  km  and  more  from  the  epicenter.  Computations  and  representations 
are  identical  as  the  ones  of  Figure  9  on  the  aftershock. 

For  this  small  foreshock  the  corner  frequency  is  of  the  order  of  7  to  9  Hz, 
to  be  compared  to  3  .  6  Hz  for  the  after  shock  (  Ml  =  4.  0  )  and  around 
2  Hz  for  the  main  shock  (  Ml  =  4.  8  ). 

19 

The  seismic  moment  is  of  the  order  of  some  10  dyne  .cm. 

But  essentialy  ,  the  spectra  definition  particularly  for  Lg}  waves  and  even 
Lg2  is  here  rather  poor  and  do  not  guarantee  a  reasonable  fiability  on  the 
re  suits . 


Conclusion  : 

This  attempt  to  evaluate  the  source  amplitude  spectrum  of  the  earthquake 
by  processing  the  Lg  waves  train  recorded  at  large  distances  seems  to  be 
reliable  as  shown  on  the  Remiremont  seismic  crisis. 

It  is  necessary  to  apply  this  method  on  more  different  cases  but  we 
notice  already  that  as  long  as  the  Lg  waves  are  well  defined,  they  bring 
information  from  the  source  which  lead  to  approach  its  main  parameters. 
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RESULTS 


Forward  calculations  have  been  completed  in  several  models  on  the  effects 
of  three-dimensional  structure  in  the  source  and  receiver  site  on  the  focussing 
and  defocussing  of  teleseismic  body  waves.  The  choice  of  models  and 
experiments  was  designed  to  separately  examine  the  effects  of  structure 
beneath  source  and  receiver  sites.  Although  the  models  that  were  examined 
differed  between  source  and  receiver  sites,  their  scale  lengths  and  intensity  of 
velocity  fluctuation  were  similar.  This  allowed  at  least  a  qualitative  test  of 
reciprocity  of  the  asymptotic  methods  used  in  the  forward  calculations.  The 
technique  used  in  the  forward  calculations  consisted  of  body  wave  synthesis  by 
dynamic  ray  tracing  and  Gaussian  beams.  In  experiments  that  investigated  the 
effects  of  receiver  structure,  a  plane  wave  incident  on  a  three-dimensional 
structure  was  expanded  into  Gaussian  beams.  In  experiments  that  investigated 
the  effects  of  source  structure,  the  propagator  matrix  of  dynamic  ray  tracing 
was  employed  to  connect  a  3-D  source  region  to  a  1-D  whole  earth  model 
(Cormier.  1986). 

Calculations  were  performed  (Nowack  and  Cormier,  19b5)  in  three- 
dimensional  models  of  the  upper  75  -  100  km.  of  the  crust  and  mantle  beneath 
the  NORSAR  array  (Thomson  and  Gubbins,  1982),  a  model  for  a  region  in 
northern  California  (Zandt,  1981),  and  a  mode'  generated  by  random 
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perturbations  to  a  1-D  velocity  structure  (McLaughlin  and  Anderson,  1985).  In 
both  the  NORSAR  and  Zandt  models,  azimuthal  variations  in  teleseismic 
amplitude  were  found  to  be  on  the  order  of  a  factor  of  2  and  variations  in  travel 
time  were  found  to  be  on  the  order  of  several  0. l’s  of  a  second.  These  models 
had  a  maximum  of  4  to  8%  velocity  fluctuation  over  scale  lengths  of  10  to  100 
km. 

A  significant  result  obtained  with  the  NORSAR  and  California  models  was 
that  the  scale-lengths  and  intensities  of  perturbations  were  such  that  the 
amplitude  variations  were  nearly  independent  of  frequency  (Figures  1  and  2). 
and  hence  adequately  predicted  by  simple  ray  theory.  This  result  has 
important  consequences  for  the  yield  estimation  of  underground  nuclear 
explosions  by  measurements  of  classical  body  wave  magnitudes,  m^,  versus 
broader  band  measurements  of  radiated  energy  in  the  time  and  frequency 
domain.  If  deep  seated,  broad  scale  length  (50km.  and  greater),  velocity 
anomalies  of  2%  or  more  are  a  common  occurrence  in  the  upper  mantle  of  the 
earth,  they  will  act  to  focus  and  defocus  body  waves  over  a  broad  frequency 
band.  The  focussing  and  defocussing  caused  by  these  broad  anomalies  will  be 
independent  of  frequency  and  will  thus  introduce  a  scatter  in  broader  band 
measures  of  radiated  energy  which  will  be  equivalent  to  that  seen  in  the  narrow 
band  measurement.  Focussing  and  defocussing  by  structure  in  either  the 
source  or  receiver  region  will  also  affect  the  coda  of  P  waves  if  a  portion  of  this 
coda  is  generated  in  either  region.  These  effects  may  help  explain  why  broader 
band  and  integrated  coda  measures  of  body  wave  energy  often  do  not  exhibit 
significantly  less  scatter  than  classical  m^  measurements. 

The  results  obtained  with  a  random  model  (Figures  3.4'  show  that  a  model 
having  a  maximum  velocity  fluctuation  as  small  as  0.8  per  cent  is  capable  of 
producing  caustics  and  multipaths  at  teleseismic  range.  The  production  of 
multipaths  strongly  depends  on  the  anisotropy  of  the  distribution  of  scale 
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lengths,  i.e.,  the  ratio  of  characteristic  vertical  to  horizontal  scale  length.  The 
multipaths  of  the  random  model,  however,  occurred  over  too  small  an  area  and 
were  too  closely  spaced  in  arrival  time  to  be  resolved  with  standard 
seismograph  systems  operating  in  the  0.01  to  4  Hz.  band. 

CONCLUSIONS 

Broad  scale  length  (50  to  100  km.),  deep  seated  structure  can  affect  short 
period  as  well  as  broad  band  and  coda  measures  of  radiated  seismic  energy.  Its 
effects,  however,  may  be  easily  correctable  if  a  3-D  model  of  the  source  region 
is  known  from  block  inversion  of  travel  time  residuals.  A  resolvable  block  size 
of  about  20  km  may  be  all  that  is  necessary  to  formulate  corrections  based  on 
azimuth  of  teleseismic  station  and  source  location  within  a  test  site. 

Calculations  with  a  random  velocity  model  demonstrated  that  a  smaller 
intensity  of  velocity  fluctuation  (0.8%)  can  produce  even  larger  amplitude 
variations  if  the  smallest  scale  length  of  fluctuation  is  on  the  order  of  10  km 
The  random  model  was  constructed  such  that  the  smallest  scale  length  of 
velocity  fluctuation  was  shorter  in  the  horizontal  direction  than  in  the  vertical 
direction.  At  teleseismic  range,  this  model  generated  isolated  caustics,  which 
were  elongated  along  the  azimuth  of  approach.  The  waveform  distortion 
associated  with  these  caustics  was  small.  The  fact  that  the  caustics  were 
generated  at  all  by  such  mild  3-D  perturbations  is  significant.  It  suggests  that 
this  type  of  synthetic  modeling  may  be  useful  in  limiting  some  of  the  attributes 
of  the  heterospectrum  of  the  earth’s  lithosphere. 

In  future  work,  the  model  of  Taylor  (1983)  for  the  Nevada  test  site  will  be 
investigated.  This  model  has  a  deep  high  velocity  anomaly  situated  to  the 
northeast  of  Pahute  Mesa,  which  may  account  for  the  deep  minimum  seen  in 
body  wave  amplitudes  of  Pahute  tests  recorded  at  teleseismic  stations  to  the 
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northeast  (Lay  et  al.,  1983).  The  patching  technique  developed  for  connecting 
3-D  models  to  1-D  radially,  symmetric  models  of  the  whole  earth  should  aid  in 
investigating  any  fine  structure  of  teleseismic  amplitudes  predicted  by  3-D 
models  of  the  source  region. 
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Figure  1  Variations  in  amplitude  for  a  vertically  incident  plane  wave  on  the 
NORSAR  model.  Results  for  several  frequencies  are  shown. 
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Figure  2  Synthetic  seismograms  constructed  by  superposition  of  Gaussian 
beams  for  stations  at  70°  and  variable  azimuths  from  an  explosive 
point  sou  ice  at  a  constant  location  SO  at  9.6  km.  depth  in  the  Zandt 
model.  Seismograms  are  shown  for  a  broad  band,  WWSSN  short  pei  iod 
and  long  period  instrument  response. 
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Figure  3  Rav  end  Doints  near  a  station  at  7T)°  distance  from  source  slO-  in 
the  random  model.  The  locations  of  a  profile  of  stations  in  the  vicinity 
of  a  caustic  intersection  is  given  by  the  tabled  points  r400-  to  r400+. 
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Synthetic  !J  waves  for  ttie  profile  of  stations  in  the  vicinity  of  the 
caustic  shown  in  the  plot  of  ray  end  points  in  Figure  3 


REGIONAL  EVENT  DETECTION  USING  THE  NORESS  ARRAY 


F.  RINGDAL 

NTNF/NORSAR,  Post  Box  51,  N-2007  Kjeller,  Norway 


INTRODUCTION 


Since  January  1985  data  from  the  small-aperture  array  NORESS  in  Norway  have 
been  processed  in  real  time  at  the  NORSAR  data  center  at  Kjeller.  The  data 
used  in  the  detection  processing  comprise  25  SPZ  channels,  deployed  over  an 
area  3  km  in  aperture  and  sampled  at  a  40  Hz  rate.  The  detection  algorithm 
has  been  described  by  Mykkeltveit  and  Bungum  (1984),  and  briefly  consists  of 

Digital  narrow-band  filtering  (six  filter) 

Beamforming  (conventional  and  incoherent) 

STA/LTA  detector  applied  to  each  beam 
Frequency-wavenumber  analysis  of  detected  signals 
Association  of  regional  phases  to  aid  in  locating  events. 

Preliminary  results  from  the  NORESS  processing  have  earlier  been  presented 
in  NORSAR  Semiannual  Technical  Summaries  (SATS).  In  this  paper,  we  discuss 
in  particular  the  automatic  detection  performance  for  events  at  regional 
distances,  and  the  spectral  characteristics  of  signal  and  noise  at  very 
high  frequencies. 


REGIONAL  DETECTABILITY 


An  assessment  of  the  NORESS  detection  capability  at  regional  distances  has 
been  obtained  by  comparing  the  automatic  NORESS  detector  output  to  the 
bulletins  produced  on  the  basis  of  local  seismic  networks  in  Fennoscandia . 
In  particular,  we  have  used  as  a  data  base  the  catalogue  of  seismic  events 
in  Northern  Europe  regularly  compiled  at  the  University  of  Helsinki. 
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The  time  period  covered  by  this  study  is  the  6-month  interval  April- 
September  1985,  during  which  the  RONAPP  processor  was  operated  with  a  fixed 
beam  deployment  (re.  NORSAR  SATS  84—85).  A  total  of  477  events  reported  in 
the  Helsinki  catalogue,  with  local  magnitudes  in  the  range  1.7-3. 3,  were 
cross-checked  with  the  NORESS  detection  list.  The  epicentral  distances  from 
NORESS  ranged  from  500-1500  km.  For  each  event,  the  expected  arrival  times 
at  NORESS  for  P,  Sn  and  Lg  were  computed,  using  standard  travel  time 
tables,  and  matched  to  NORESS  detection  entries. 

Fig.  1  shows  the  magnitude-distance  distribution  of  the  events  in  the  data 
base.  Events  detected  at  NORESS  are  shown  as  crosses,  non-detections  are 
indicated  as  triangles.  In  this  figure,  "detection"  means  that  at  least  one 
phase  (P,  Sn  or  Lg)  was  reported  by  NORESS. 

A  summary  of  the  statistics  on  automatic  detection  of  primary  vs.  secondary 
phases  is  given  in  Table  1.  We  note  that,  at  low  magnitudes,  many  events 
are  detected  only  on  secondary  phases.  It  is  also  noteworthy  that  several 
events,  even  in  the  higher  magnitude  range,  are  detected  as  P-phases  only. 
However,  visual  inspection  of  the  signal  traces  shows  that  in  virtually  all 
of  these  cases  an  Lg  phase  may  be  identified  by  the  analyst.  Thus,  the  lack 
of  secondary  phase  detections  is  a  problem  within  the  automatic  processor 
that  requires  improvements  in  the  algorithms  in  order  to  extract  emergent 
phases  in  the  coda  of  a  preceding  P-phase. 

Table  1 

Summary  of  automatic  NORESS  detection  statistics  for  the  regional  data  base. 


Mr. 

1. 5-2.0 

2. 1-2.5 

2. 6-3.0 

3. 1-3.5 

P  detection  only 

13 

39 

16 

0 

P  +  secondary  phase* 

16 

105 

30 

3 

Secondary  phase  only* 

28 

88 

2 

0 

No  detection 

48 

87 

2 

0 

Total 

105 

319 

50 

3 

*  "Secondary  phase"  meaning  Sn  or  Lg  (or  both)  detected. 
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Fig.  1  shows  that  there  is  only  a  slight  degradation  in  detection  capabi¬ 
lity  with  increasing  distance  for  the  range  considered.  As  an  initial  esti¬ 
mate  of  detection  thresholds,  we  have  therefore  combined  the  data  in  the 
distance  range  700-1400  km  and  estimated  detection  thresholds  as  shown  in 
Fig.  2  (detection  on  either  P,  Sn  or  Lg)  and  Fig.  3  (P-detection  only).  The 
method  described  by  Ringdal  (1975)  has  been  applied. 

From  Fig.  2  we  observe  that  the  50  and  90  per  cent  regional  detection 
thresholds  are  close  to  Ml  =  1.9  and  2.5,  respectively.  When  only  P  phase 
detections  are  counted  (Fig.  3),  the  respective  thresholds  are  Ml  =  2.3  and 
2.7. 

It  would  clearly  be  desirable  to  tie  these  thresholds  to  the  global  m^ 
scale.  It  has,  however,  not  been  possible  to  do  this  for  the  present  data 
set,  since  all  of  the  477  reference  events  are  much  too  small  to  have  any 
teleseismic  detection.  Nevertheless,  the  local  magnitude  scale  in  question 
has  been  developed  so  as  to  be  consistent  with  world-wide  m^,,  and  the  dif¬ 
ferences  are  not  thought  to  be  significant.  This  topic  will  be  subject  to 
further  study. 

As  a  final  note,  we  remark  that  the  large  majority  of  reference  events  are 
mining  explosions,  mostly  from  mines  in  Western  Russia.  We  have  not  yet 
attempted  systematically  to  compare  the  detectability  of  these  explosions 
to  that  of  the  (very  few)  earthquakes  in  the  data  base,  but  initial  studies 
do  not  indicate  major  differences  for  the  two  source  types. 

HIGH-FREQUENCY  STUDIES 

The  recently  installed  high-frequency  recording  system  (HFSE)  at  NORESS  has 
provided  a  unique  opportunity  to  study  noise  and  signal  characteristics  at 
frequencies  up  to  at  least  50  Hz.  These  studies  have  only  begun,  but  it  is 
alredy  apparent  that  much  important  information  can  be  extracted  from  the 
high  frequency  part  of  the  spectrum.  Based  on  about  50  regional  events  ana¬ 
lyzed  so  far,  it  appears  that  the  signal-to-noise  ratio  of  P  and  Sn  phases 
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does  not  show  any  appreciable  decay  with  increasing  frequency  up  to  at 
least  50  Hz  for  events  within  about  500  km  epicentral  distance.  At  greater 
distances,  the  signal  spectra  drop  more  rapidly  than  the  noise  spectra,  but 
even  at  1000  km  distance,  the  P-wave  of  Ml  =  3.0  events  have  significant 
signal-to-noise  ratio  up  to  about  25-30  Hz. 

An  example  of  a  high-frequency  recording  is  shown  in  Fig.  4,  corresponding 
to  an  Ml  =  5.0  earthquake  off  the  west  coast  of  Norway  on  February  6,  1986 
(distance  =  417  km).  The  unfiltered  record  shows  the  expected  amplitude 
pattern,  i.e.,  Lg  as  the  dominant  phase,  Pg  much  larger  than  Pn.  The  pic¬ 
ture  changes  dramatically  when  considering  the  high-frequency  part  of  the 
record.  In  the  filter  band  30-50  Hz,  the  Pn  and  Sn  phases  dominate  the 
seismogram,  and  the  Pg  and  Lg  phases  are  not  even  visible. 

Fig.  5  shows  HFSE  spectra  from  the  vertical  component  for  the  same  event.  We 
see  that  Lg  exceeds  the  preceding  noise  (which  in  fact  is  the  Sn  coda)  only 
up  to  about  10  Hz,  whereas  Pn  has  large  SNR  over  the  entire  frequency  band. 

In  conclusion,  the  NORESS  regional  detection  capability  appears  to  be 
Ml  =  2.5  or  better  out  to  at  least  1500  km.  Many  small  events  are  detected 
only  on  secondary  phases,  and  the  further  improvement  of  automatic  detec¬ 
tion  of  such  phases  is  important.  At  distances  up  to  500  km,  considerable 
improvements  in  detection  capability  are  possible  by  taking  advantage  of 
the  high  frequencies  which  propagate  very  efficiently  in  this  distance 
range.  The  high  frequency  band  is  also  potentially  valuable  for  improved 
phase  identification,  especially  to  separate  Pn  from  Pg  and  Sn  from  Lg. 
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NORESS  REGIONAL  DETECTION 

AZIMUTH  RANGE  15-175  DEGREES 


Distribution  of  489  regional  events  in  the  data  base  as  a  function 
of  epicentral  distance  from  NORESS  and  local  magnitude  Ml.  Crosses 
denote  events  detected  automatically  at  NORESS  (either  P,  Sn  or  Lg 
phase),  whereas  non-detected  events  are  marked  as  triangles. 
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Time  domain  plots  of  the  HFSE  recordings  on  the  SPZ  channel  at 
NORESS  for  an  Ml  *  5.0  earthquake  off  the  west  coast  of  Norway 
(distance  -  417  km).  The  upper  part  covers  the  entire  wavetrain 
(unfiltered  and  in  two  filter  bands  as  indicated).  The  bottom  part 
is  an  expanded  view  of  the  P  window.  Note  the  prominence  of  Pn  and 
Sn  in  the  high-frequency  band. 
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.  5  Spectral  plot  of  the  Pn  and  Lg  phases  for  the 
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SUMMARY 

Six  events  showing  teieseismic  P  coda  recorded  at  the  NORESS  array  in  the  distance 
range  35°-85°have  been  analysed.  The  events  were  two  Semipalatinsk  presumed  explosions, 
two  deep  focus  400  km)  earthquakes  in  the  Japan-Bonin  trench,  and  two  crustal  earth¬ 
quakes  with  focal  depth  <  15  km  in  Alaska  and  Yunnan,  China.  The  analysis  methods  were 
measurement  of  coda  power  relative  to  the  energy  in  first  P,  determination  of  the  decay  rate  in 
time  of  the  coda  power  spectrum,  frequency-wavenumber  analysis,  and  measurement  of  the 
wavenumber  spectrum.  The  model  tested  states  that  the  coda  consists  of  two  components 
produced  by  scattering  of  Lg  to  P  in  the  crust  near  the  source  and  P  to  Lg  near  the  receiver. 
The  near  source  component  will  not  be  present  for  deep  focus  events.  The  model  predicts:  ( 1 ) 
crustal  earthquakes  have  highest  power  in  the  coda  relative  to  first  P  and  deep  focus  events 
lowest;  (2)  coda  decays  exponentially;  (3)  the  near  source  component  has  the  same  phase 
velocity  and  azimuth  as  first  P,  the  near  receiver  component  has  the  phase  velocity  of  Lg  in  the 
receiver  region  and  random  azimuth;  (4)  the  differences  predicted  in  (1 )  are  due  to  differences 
in  the  strength  of  the  near  source  component.  The  tests  (1)  and  (2)  are  considered  weak 
tests  of  the  model,  but  (3)  and  (4)  seem  to  be  quite  demanding.  The  model  passes  all  four 
tests  on  the  events  examined. 
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INTRODUCTION 


Measurements  on  teleseismic  and  regional  coda  have  been  proposed  for  the  estimation  of 
magnitude  (e.g.  Bullitt  and  Cormier,  1984)  and  as  a  discriminant  (Dainty  and  Harris,  1985).  The 
nature  of  the  coda  is  not  well  understood,  however;  it  apparently  consists  of  scattered  waves 
(Aki  and  Chouet,  1 975)  but  the  wave  types  involved  and  whether  single  or  multiple  scattering 
is  important  are  subjects  of  controversy  (Dainty,  1 985).  This  study  is  an  attempt  to  determine 
the  wave  types  present  in  coda  and  the  location  of  the  scatterers  for  events  recorded  at  the 
NORESS  regional  array.  To  date,  the  work  has  encompasssed  six  teleseismic  events  (Table  1); 
coda  from  20  to  1 20  seconds  after  first  P  has  been  examined.  The  distance  to  these  events 
ranges  from  35°  to  85°.  The  events  consist  of  two  presumed  nuclear  explosions  from  the 
Semipaiatinsk  test  site,  two  deep  focus  earthquakes  from  the  Japan-Bonin  trench  and  two 
shallow  focus  (crustal)  earthquakes. 


Table  1 


Date 

OT 

(UTC) 

Latitude 

Longitude 

Depth 

(km) 

mb 

Agency 

Location 

1984  Nov  17 

18:27:13.1 

28.28694 

1 39.8499: 

465 

5.3 

PDE 

Bonin 

1985  Feb  10 

3:27:7.6 

49.87794 

78.81 69: 

0 

5.9 

EDR 

Semipaiatinsk 

1985  Mar  9 

14:8:4.1 

66.21 594 

150.0639^ 

12 

5.9 

PDE 

Alaska 

1985  Apr  10 

16:26:18.6 

29.97994 

138.7909 

398 

5.8 

PDE 

Honsh!' 

1985  Apr  18 

5:52:5 2.7 

25.89894 

102.8709 

5 

5.7 

PDE 

Yunnan 

1985  Apr  25 

0:57:6.5 

49.907*94 

78.9329 

0 

5.9 

PDE 

Semipaiatinsk 

ANALYSIS  AND  MODEL 


The  principal  methods  of  analysis  are  determination  of  coda  and  first  P  power  as  a  func¬ 
tion  of  frequency  and  time,  standard  frequency-wavenumber  analysis  and  wavenumber  spectra. 
A  brief  description  of  these  methods  is  given  in  Dainty  and  Harris  (1985);  this  discussion  will 
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emphasise  additions  to  their  description.  Power  as  a  function  of  time  and  frequency  is  found  by 
taking  the  power  spectrum  in  overlapping  five  second  windows  and  averaging  over  a  one 
octave  frequency  band.  Frequency-wavenumber  spectra  for  linear  wavenumber  North  and  East 
are  found  for  five  second  windows  by  standard  techniques  at  a  specified  frequency,  with 
averaging  in  the  coda.  The  averaging  is  accomplished  by  determining  frequency-wavenumber 
spectra  in  contiguous  windows  over  a  specified  time  interval,  correcting  for  the  coda  decay 
found  from  analysis  of  coda  power  as  a  function  of  frequency  and  time,  and  averaging  the 
corrected  spectra.  The  wavenumber  spectrum  estimates  the  power  in  the  wavefield  as  a  func¬ 
tion  of  linear  wavenumber  without  regard  to  azimuth  at  a  specified  frequency  f:  it  is  equivalent 
to  Integrating  a  frequency-wavenumber  spectrum  around  a  circle  centered  on  the  wavenumber 
origin  (Dainty  and  Harris,  1 985).  Averaging  in  the  coda  is  accomplished  by  the  same  technique 
used  for  frequency-wavenumber  spectra.  Since  the  linear  wavenumber  k  =  f/c,  where  c  is  the 
phase  velocity,  the  wavenumber  spectrum  has  the  potential  of  identifying  components  of  the 
wavefield  by  their  phase  velocity. 

The  results  of  these  analyses  have  been  compared  with  a  model  of  teleseismic  P  coda.  In 
this  model  the  coda  consists  of  two  dominant  components  due  respectively  to  scattering  of  Lg 
to  teleseismic  P  near  the  source  and  scattering  of  teleseismic  P  to  Lg  near  the  receiver.  If  sin¬ 
gle  scattering  is  assumed,  the  coda  decays  exponentially  (Dainty  and  Harris,  1985).  The 
decay  may  be  characterised  by  a  Q,  and  the  ratio  of  the  coda  power  to  the  energy  in  first  P  is 
proportional  to  the  side  scattering  turbidity  (side  scattering  cross-section  per  unit  volume)  in 
the  source  and  receiver  regions.  This  implies  that  the  power  ratio  will  depend  on  the  source 
type  and  depth:  deep  focus  earthquakes  will  have  a  small  ratio  (weak  coda)  due  to  lack  of 
scattering  near  the  source,  whereas  crustal  earthquakes  will  have  strong  coda  due  to  the  high 
level  of  Lg  avilable  for  scattering  near  the  source.  Explosions  will  have  a  contribution  from  the 
source  region,  but  not  as  strong  as  crustal  earthquakes. 

Frequency-wavenumber  analysis  and  wavenumber  spectra  should  allow  a  quantitative 
identification  of  the  two  components  of  the  coda.  For  explosions  and  crustal  earthquakes 
frequency-wavenumber  spectra  of  the  coda  should  show  a  peak  at  the  same  phase  velocity 
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and  azimuth  as  first  P,  representing  energy  scattered  near  the  source  and  travelling  to  the 
receiver  as  P.  Energy  scattered  near  the  receiver  should  have  random  azimuth  and  hence  might 
not  be  evident  on  frequency-wavenumber  analysis,  but  wavenumber  spectra  should  show  a 
peak  at  high  wavenumber  (low  phase  velocity,  c  ~  3.5  km/s)  corresponding  to  energy  travelling 
as  Lg.  Explosion  and  crustal  earthquake  wavenumber  spectra  should  also  show  a  peak  at  low 
wavenumber  corresponding  to  the  high  phase  velocity  component  from  the  source  region;  the 
relative  size  of  the  peaks  will  indicate  the  relative  contribution  of  scattering  near  the  source 
and  near  the  receiver. 


RESULTS  AND  DISCUSSION 

The  analyses  discussed  above  were  carried  out  on  the  events  of  Table  1  for  the  fre¬ 
quency  range  0.9  -  7.2  Hz.  Measurements  of  coda  decay  and  coda  power  to  first  P  ratio  were 
carried  out  for  frequencies  of  0.9,  1.8,  3.6  and  7.2  Hz,  provided  signal  to  noise  ratio  was  ade¬ 
quate,  for  the  Bonin,  Semipalatinsk  (Feb  1 0),  Alaska  and  Yunnan  events.  Q  increased  with  fre¬ 
quency,  ranging  from  160  at  0.9  Hz  for  Yunnan  to  220<"<  at  7.2  Hz  for  Bonin.  If  the  decay 
results  are  cast  in  terms  of  the  apparent  total  turbidity  (Dainty,  1985),  this  quantity  is  approx¬ 
imately  independent  of  frequency  and  has  an  average  value  of  0.005  km-1  (mean  free  path  of 
»200  km).  The  Yunnan  event  has  an  apparent  total  turbidity  of  0.01  km-1,  presumably  reflect¬ 
ing  greater  attenuation  in  the  source  region.  The  ratio  of  coda  power  to  first  P  energy  shows  a 
stronger  coda  relative  to  first  P  for  the  two  crustal  earthquakes,  as  expected,  but  is  rather 
variable,  probably  because  of  focussing-defocussing  effects  on  first  P.  These  results  are  con¬ 
sistent  with  the  model,  but  could  likely  be  consistent  with  other  models  too. 

Results  from  frequency-wavenumber  analysis  and  wavenumber  spectra  were  more  defina- 
tive.  In  both  methods  the  analyses  were  carried  out  at  the  peak  power  frequency,  which 
ranged  between  1  and  2  Hz.  Frequency-wavenumber  analysis  for  the  explosions  and  crustal 
earthquakes  shows  a  strong  component  of  the  coda  with  phase  velocity  and  azimuth  similar  to 
first  P;  the  deep  focus  earthquakes  had  a  component  with  a  similar  azimuth  to  first  P,  but  much 
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slower  (5-6  km/s)  phase  velocity.  These  results  are  in  accord  with  the  model.  The 
wavenumber  spectra  provide  even  stronger  evidence  for  the  model.  Figures  1-3  show  such 
spectra  for  the  Bonin  (deep  focus),  Semipalatinsk  (Feb  10;  explosion)  and  Yunnan  (crustal) 
events,  for  both  first  P  and  coda.  There  is  a  biasing  effect  towards  high  wavenumber  (low 
phase  velocity)  due  to  the  finite  aperture  (3  km)  of  the  NORESS  array,  which  may  be  assessed 
by  examining  the  first  P  peak:  true  phase  velocities  for  first  P  are  between  12  and  16  km/s. 
The  effect  should  be  less  for  lower  phase  velocities.  The  deep  focus  event  (Figure  1 )  shows  a 
pronounced  shift  from  high  velocity  (low  wavenumber)  energy  in  first  P  to  low  phase  velocities 
typical  of  Lg  in  the  coda,  indicating  that  the  coda  consists  of  energy  scattered  near  the 
receiver.  The  explosion  (Figure  2)  however  retains  a  component  of  similar  phase  velocity  to 
first  P  in  the  coda,  of  about  equal  amplitude  to  the  low  velocity  Lg  component.  The  crustal 
event  (Figure  3)  has  a  coda  dominated  by  high  phase  velocity  energy  that  must  have  been 
scattered  near  the  source.  Lg  to  P  scattering  near  the  source  is  attractive  because  it  would 
explain  the  difference  between  the  explosion  (generates  low  Lg)  and  the  crustal  earthquake 
(generates  high  Lg). 
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BONIN  DEEP  FOCUS  1  2  HZ  15-20  SEC 


BONIN  DEEP  FOCUS  1.2  HZ  40  -  120  SEC 


Figure  1 .  Wavenumber  spectra  for  the  Bonin  event.  First  P  is  shown  at  the  top  with  coda  from 
25  to  105  sec  after  first  P  at  the  bottom.  Frequency  is  shown  on  the  Figure; 
selected  phase  velocities  are  indicated. 
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S  EM  IP  ALA  TINSK  1.8  HZ  15-20  SEC 


SEMIPALATINSK  1 .8  HZ  40  -  120  SEC 


Figure  2.  As  for  Figure  1,  but  for  Semipalatinsk  event  of  Feb  10. 
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YUNNAN  1.1  HZ  20  -  25  SEC 


„  YUNNAN  1.1  HZ  40  -  120  SEC 


Figure  3.  As  for  Figure  1,  but  for  Yunnan  event.  Coda  wavenumoer  spectrum  is  for  20  to  100 
sec  after  first  P. 
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PROPAGATION  CHARACTERISTICS  OF  REGIONAL  PHASES 


AND  NORESS  REAL  TIME  PROCESSING  PERFORMANCE 


S.  Mykkeltveit  and  T.  Kvaarna 
NTNF/NORSAR,  Post  box  51,  N-2007  Kjeller,  Norway 


The  performance  of  the  RONAPP  real  time  processing  package  (Mykkeltveit  and 
Bungura,  1984)  for  regional  events  recorded  on  the  NORESS  array  is  currently 
being  evaluated.  Detailed  knowledge  of  propagation  path  effects  must  be 
obtained  for  the  purpose  of  calibrating  the  real  time  processor.  Thie 
contribution  gives  i)  the  results  from  a  general  mapping  study  aimed  at 
obtaining  such  calibration  data  and  ii)  the  results  from  an  investigation 
of  RONAPP  performance  for  events  in  Western  Norway  and  the  North  Sea. 

PROPAGATION  CHARACTERISTICS  OF  REGIONAL  PHASES  RECORDED  AT  NORSAR 


Previous  investigations  (Kennett  and  Mykkeltveit,  1984;  Kennett  et  al, 

1985)  limited  to  specific  regions  have  shown  that  prominent  geological 
features  exist  within  regional  distance  from  NORESS  that  strongly  influence 
the  characteristics  of  secondary  phases.  We  have  recently  completed  a  study 
utilizing  data  from  the  NORSAR  data  archive,  corresponding  to  regional 
events  at  all  azimuths.  For  the  170  events  in  this  data  base,  the  seismo¬ 
grams  were  band-pass  filtered  in  three  different  frequency  bands,  centered 
on  1.0,  2.0  and  4.0  Hz.  For  each  of  these  bands,  arrivals  corresponding  to 
distinct  phase  onsets  were  picked  by  an  STA/LTA  detector.  In  this  manner, 
our  off-line  analysis  ties  in  with  the  on-line  processing  of  NORESS  data 
with  respect  to  detection  of  the  regional  phases. 

The  primary  aim  of  the  analysis  was  to  characterize  the  relative  strength 
of  the  Sn  and  Lg  phases.  The  case  of  dominant  Lg  (over  Sn)  was  subdivided 
into  three  different  codes  as  follows: 

(1)  Only  Lg  can  be  seen  in  the  S  wavetrain. 

(2)  Lg  is  the  dominating  phase,  but  Sn  can  be  discerned. 

(3)  Lg  is  still  the  dominating  phase,  but  Sn  is  not  very  much  smaller. 
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Ail  equivalent  subdivision,  also  on  a  three-point  scale,  was  adopted  for  the 
case  of  dominant  Sn.  Finally,  a  code  was  assigned  to  the  case  of  comparable 
Lg  and  Sn  amplitudes. 

Fig.  1  shows  the  results  derived  from  analysis  of  the  data  in  our  data 
base.  In  this  figure,  we  assign  the  character  of  the  S  wavetrain  as 
observed  in  the  NORSAR  seismometer  record  to  the  source  area.  We  have  con¬ 
fined  ourselves  to  the  use  of  three  codes  only  in  this  figure,  but  color- 
coded  maps  including  all  7  codes  defined  above  have  also  been  prepared. 
These  maps  offer  more  details  and  often  show  gradual  transitions  from  domi¬ 
nating  Lg  to  dominating  Sn  in  certain  directions.  Inspection  of  Fig.  1  and 
a  tectonic  map  of  the  same  area  leads  to  the  following  conclusions: 

Lg  wave  propagation  is  very  efficient  for  propagation  paths  that  are 
confined  to  continental  shield  type  structures  (dominant  to  the  east 
of  NORSAR)  not  intersected  by  major  tectonic  units.  In  these 
directions,  1  Hz  Lg  waves  propagate  to  20°  from  the  source,  while  at 
4  Hz  Lg  is  seen  out  to  about  10°  from  NORSAR. 

No  Lg  is  observed,  or  Sn  dominates  over  Lg  for  propagation  paths 
that  cross  significant  tectonic  features  like  graben  structures  (in 
the  North  Sea)  and  continental  margins. 

-  Particularly  at  4  Hz  there  is  a  close  connection  between  some  of  the 

finer  tectonic  details  and  the  propagation  characteristics.  Since  Lg 
is  a  phase  mostly  composed  of  waves  trapped  in  the  crust,  this  phase 
is  likely  to  be  highly  susceptible  to  changes  in  the  crustal 
structure . 

RONAPP  PERFORMANCE  FOR  EVENTS  IN  WESTERN  NORWAY  AND  THE  NORTH  SEA 

For  an  evaluation  of  RONAPP  we  would  like  to  find  a  region  with  events  that 
have  been  reliably  located  by  another  agency,  based  on  a  local  network.  One 
such  region  within  regional  range  from  NORESS  is  Western  Norway  with  the 
adjacent  North  Sea.  During  the  period  April-September  1985  a  total  of  230 
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events  in  this  region  were  reported  by  the  University  of  Bergen.  These 
events  are  typically  earthquakes  offshore  and  explosions  onshore.  The 
magnitudes  are  in  the  range  Ml  1.5-3. 5,  and  the  epicentral  distances  rela¬ 
tive  to  NORESS  range  from  200  to  nearly  700  km. 

The  seismograms  from  these  events  are  rich  in  regional  phases,  and  quite 
frequently  Pn,  Pg,  Sn  and  Lg  are  clearly  seen  in  the  records.  An  example 
is  shown  in  Fig.  2.  As  can  be  seen  from  Fig.  1,  Lg  is  dominant  for  onshore 
locations  in  Western  Norway,  whereas  the  importance  of  Sn  increases  with 
increasing  frequency  for  offshore  locations.  For  higher  frequencies  than 
those  in  Fig.  1,  the  importance  of  Sn  increases  further.  Data  from  a  high- 
frequency  system  (125  Hz  sampling  rate)  recently  installed  within  NORESS 
shows  considerable  Sn  energy  in  the  band  30-50  Hz  from  an  earthquake  (PDE 
mb  5.0)  in  this  region  at  a  distance  of  417  km  from  NORESS  (see  Ringdal, 
this  volume). 

The  real  time  location  by  RONAPP  of  a  regional  event  results  from  con¬ 
secutive  detections  of  P-  and  S-type  arrivals  with  a  common  azimuth.  The 
range  is  estimated  from  the  travel-time  difference  between  the  first  P 
phase  (assumed  to  be  Pn)  and  the  strongest  phase  in  the  S  wavetrain,  which 
is  assumed  to  be  Lg,  propagating  at  a  constant  group  velocity  of  3.5  km/s. 
Fig.  3  shows  the  differences  between  P  and  Lg  azimuths  for  the  class  of 
events  that  was  also  automatically  located  by  RONAPP,  based  on  the  detec¬ 
tion  of  both  a  P-  and  S-type  arrival.  The  median  of  the  absolute  values  of 
these  differences  is  5  degrees.  This  is  an  estimate  then  of  the  consistency 
that  should  be  expected  among  the  azimuths  automatically  determined  for  the 
various  phases.  Fig.  4  shows  the  differences  between  the  ranges  estimated 
by  the  NORESS  online  processing  and  the  ones  corresponding  to  the  network 
solutions.  As  can  be  seen  from  this  figure,  there  is  a  tendency  by  the 
NORESS  event  location  procedure  to  underestimate  the  range  for  events  in 
this  region.  The  general  tendency  of  underestimation  of  range  can  be 
explained  by  i)  the  fact  that  the  travel-time  model  used  in  RONAPP  is  based 
on  a  40  km  thick  crust,  whereas  the  crust  is  typically  30  km  thick  in 
Western  Norway  and  ii)  that  for  approximately  30  events  we  found  reason  to 
believe  that  the  Pn  phase  had  gone  undetected,  and  the  P  phase  detected  is 
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most  likely  to  be  Pg.  Fig.  4  shows  that  there  are  approximately  20  events 
for  which  the  range  difference  exceeds  40  km.  By  inspecting  the  relevant 
plots,  it  was  found  that  the  cases  with  large  underestimation  of  the  range 
are  related  to  mistaking  Sn  for  Lg.  With  proper  phase  assignments,  the 
ranges  for  thece  9  events  are  correctly  determined.  The  main  reason  for 
overestimating  the  range  is  that  the  strongest  detection  n  the  S  wavetrain 
is  in  the  Lg  coda. 

Fig.  5  shows  deviations  from  expected  arrival  times  (based  on  the  network 
locations)  for  all  detections  by  RONAPP  that  can  be  associated  with  the  230 
events  (the  criterion  is  that  the  azimuth  estimated  for  the  detected  phase 
must  be  within  30°  of  the  azimuth  from  the  network  location).  These  detec¬ 
tions  are  classified  as  P  detections  if  the  estimated  phase  velocity 
exceeds  6  km/s,  and  as  either  Sn  or  Lg  (depending  on  arrival  time  relative 
to  expected  arrival  times  for  Sn  and  Lg)  if  the  phase  velocity  is  below 
this  value.  As  is  indicated  by  this  figure,  the  class  of  P  detections  is 
made  up  of  slightly  early  Pn  detections  (can  be  corrected  by  travel  time 
model  modifications)  and  Pg  arrivals  typically  4  or  3  seconds  after  Pn.  The 
Lg  onset  detections  are  typically  early  by  about  3  seconds,  and  there  is  a 
fair  amount  of  Lg  coda  detections.  These  results  suggest  that  the  group 
velocity  of  3.50  km/s  used  by  RONAPP  for  the  Lg  onset  in  the  range  estima¬ 
tion  should  be  changed  to  approximately  3.60  km/s  for  this  region. 
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1  Characterization  of  Lg  versus  Sn  propagation  efficiency  for  the 
three  different  filter  bands.  The  shading  codes  are  as  follows:  (1) 
Lg  dominant  phase  in  the  S  wavetrain;  (2)  Lg  and  Sn  comparable  and 
(3)  Sn  dominant.  Note  that  the  codes  are  assigned  to  the  source 
area,  and  all  characterizations  are  related  to  seismograms  as  they 
appear  on  NORSAR  (filled  circle)  seismometers. 


Fig.  2  Individual  NORESS  traces  for  one  of  the  events  (M^,  2.1,  distance 
350  km)  reported  by  the  University  of  Bergen.  The  panel  covers  61 
minutes  of  bandpass  filtered  records  (3-5  Hz).  The  four  arrows  on 
top  of  the  records  indicate  detections  by  the  RONAPP  real  time  pro¬ 
cessing  package.  The  broad  arrows  correspond  to  Pn  and  Lg  arrivals, 
which  have  been  combined  to  locate  the  evnet  at  the  geographical 
position  given  on  top  of  the  figure.  The  two  additional  detections 
correspond  to  Pg  and  Sn  arrivals.  P  and  Lg  beams  are  shown  in  the 
bottom  part  of  the  figure. 
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NORESS  EVALUATION 

NORESS  P-LC  AZIMUTH  DIFFERENCES 
DATA  BASE  230  WESTERS  NORWAY  EVENTS 


AZIMUTH  DIFFERENCE  (DEC) 


Fig.  3  The  histogram  shows  the  differences  between  the  azimuths  automati¬ 
cally  estimated  for  the  P  and  Lg  phases  by  the  routine,  real  time 
processing  of  NORESS  data,  for  the  events  reported  by  the  University 
of  Bergen. 


NORESS  EVALUATION 

RANCE  DIFFERENCES  (  NORESS  MINUS  NETWORK  ) 
DATA  BASE  230  WESTERN  NORWAY  EVENTS 


RANCE  DIFFERENCE  (KM) 


Fig.  4  The  histogram  shows  the  differences  between  the  ranges  estimated  by 
the  NORESS  online  processing  and  the  ones  corresponding  to  the  net¬ 
work  solutions.  35  events  with  known  location  at  an  explosion  site 
at  range  290  km,  for  which  this  difference  was  between  -10  and 
-20  km,  are  not  included  in  the  figure. 
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NORESS  P  DETECTIONS 
DATA  BASE  230  WESTERN  NORWAY  EVENTS 


DEVIATION  FROM  EXPECTED  ARRIVAL  TIME  (SEC) 


NORESS  SN  DETECTIONS 
DATA  BASE  230  WESTERN  NORWAY  EVENTS 


NORESS  LG  DETECTIONS 
DATA  BASE  230  WESTERN  NORWAY  EVENTS 


DEVIATION  FROM  EXPECTED  ARRIVAL  TIME  (SEC) 


Fig.  5  Deviations  from  expected  arrival  times  for  P  (top),  Sn  (middle)  and 
Lg  (bottom)  for  NORESS  detections  corresponding  to  the  230  Western 
Norway  events. 
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WAVEFIELD  DECOMPOSITION  USING  ML-PROBABILITIES 


IN  MODELLING  SINGLE-SITE  3-COMPONENT  RECORDS:  PRACTICAL  APPLICATIONS 
ON  LOCATIONS  OF  EVENTS  AND  SCATTERING 


E.S.  Husebye^,  A.  Christoff ersson^) ,  S.F.  Ingate^)  and  B.O.  Ruud*) 
1)  NTNF/NORSAR,  Kjeller,  Norway;  2)  univ.  of  Uppsala,  Sweden; 

3)  MIT,  ERL,  Cambridge,  USA 


Abstract.  In  this  paper  a  novel  approach  to  3-component  seismogram  analy¬ 
sis  is  presented.  The  P,  S,  Love  and  Rayleigh  model  representations  are 
similar  to  those  used  previously,  but  we  use  a  Maximum  Likelihood  (ML) 
estimator  and  rather  than,  e.g.,  Principal  Component  (PC)  estimator.  The 

reason  is  that  the  ML-approach  permits  a  model  test  in  terms  of 
2 

X  -probabilities  on  whether  a  specific  phase  or  wavelet  type  is  present  in 
the  records  as  a  function  of  time  and  azimuth.  In  addition,  apparent  angle 
of  incidence  or  apparent  velocity  is  estimated  jointly  with  probability. 

For  a  single  3-component  time-domain  seismogram,  the  covariance  matrix  pro¬ 
vides  only  6  independent  observations,  thus  restricting  analysis  to  rather 
simple  wave  models.  Exploiting  wavefield  structure  in  this  manner  appears 
very  beneficial  for  seismogram  analysis;  practical  applications  hereof  will 
be  given  in  terms  of  epicenter  determinations  at  any  distance,  locations  of 
scattering  sources  and  noise  directionality. 

1 .  INTRODUCTION 


Polarization  properties  of  scalar  and  vector  fields  have  been  topics  of 
extensive  research  in  areas  ranging  from  electromagnetics  (optics  and  radio 
transmission)  to  physical  oceanography,  seismology  and  not  at  least  seismic 
prospecting  including  VSP-prof iling .  Seismologists  have  for  many  years 
attempted  to  exploit  the  information  potential  of  three-component  records 
for  wave  propagation  modelling,  retrieval  of  structural  information  and 
source  parameters.  Such  efforts  appear  to  have  been  frustrated,  primarily 
due  to  the  technique  of  analysis  employed  -  in  other  words,  a  certain 
failure  to  produce  results  in  a  useful  manner,  at  least  for  frequencies 
above  0.2  Hz.  Therefore,  3-component  seismogram  analysis  techniques  have 
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not  become  a  viable  alternative  to  those  often  used  for  arrays  of  vertical- 
component  seismometers  such  as  frequency-wavenumber  (f-k)  analysis.  A 
disadvantage  of  arrays  is  that  only  amplitude  and  timing  information  is 
utilized,  excluding  wavefield  polarization  properties.  Our  novel  technique 
for  extracting  wavefield  structure  is  given  in  section  2,  and  then  asso¬ 
ciated  practical  analysis  results  follow. 

2.  THEORY 


The  method  for  analysis  of  3-component  data  is  based  on  the  following 
model : 


y(t)  =  A  z(t)  +  e(t)  (1) 

where  y(t)  =  [yi(t) ,y2(t)Y3(t)  ]*  is  a  vector  containing  the  observed 
3-component  recordings  at  time  t,  and  *  denotes  vector  transpose.  Simi¬ 
larly,  z(t)  is  a  k-dimensional  representation  of  the  signal  at  time  t; 
k  =  1  in  cases  of  P,  SV  and  SH  motion  and  the  A  is  a  matrix  of  unknown 
constants : 


I  Ml  •••  Mk  j 
Ml  •••  Mk  [ 

Ml  •••  Mk  ] 

relating  the  signal  to  the  observed  data.  e(t)  is  the  residual  noise  vec¬ 
tor.  Four  basic  assumptions  are  made:  (i)  z(t)  and  e(t)  are  uncorrelated; 
(ii)  e(t)  has  zero  expectation,  i.e.,  E(e(t))  =  0;  (iii)  the  components  of 
z(t)  are  linearly  independent,  i.e.,  the  signal  is  k-dimensional,  with  no 
interference  between  P,  SV  and  SH;  and  (iv)  All  moments  up  to  at  least 
second  order  exist.  These  basic  assumptions  imply  that  the  zero-lag  second 
order  moments  of  y(t)  can  be  written  as 

E(t)  =  A  <$(t)A*  +  <|»(t)  (2) 

where  £(t)  -  E(y(t)y(t)*) ;  <S(t)  =■  E(z( t  )z( t )* ) ;  and  <|»(t)  =  E(e(t)e(t)*). 
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Averaging  eq.  (2)  over  a  time  widow  Tj<t<T2  yields  the  following  structure: 
£  =  A  $  A*  +  4,  (3) 

t2  t2  t2 

where  £  =  /  £(t)  dt  ;  <$  =  /  <£(t)  dt  ;  =  j  <Kt)  dt 

Tl  T1  T1 

is  the  expected  3x3  symmetric  covariance  matrix  of  the  recordings, 
the  signal  and  noise  covariance  matrices,  respectively.  The  structure  of 
the  A  and  $  matrices  depends  on  the  particle  motion  of  the  incoming  wave 
beirg  P,  SV,  SH,  L  (Love)  or  R  (Rayleigh).  Further  details  on  the  modelling 
are  given  by  Christoff ersson  et  al,  1986  (submitted  JGRAS). 

2.1  Estimation 

The  unknown  parameters  in  equation  (3)  can  be  estimated  on  the  basis  of  the 
observed  second-order  zero-lag  moments  over  a  time  window 

1  T2 

s  =  -  •  l  y(t)  y(t)*  (4) 

T2-Ti+1  t-Ti 

There  are  several  possible  estimators  that  could  be  used.  Most  are 
based  on  some  fitting-function,  i.e.,  minimizing  some  weighted  func¬ 
tion  F  of  the  difference  between  observed  S  and  the  theoretical  E. 

Our  preference  so  far  is  for  the  Maximum  Likelihood  (ML)  estimator,  which 
is  based  on  Gaussian  assumptions  and  minimizes 

F  =  log  | £|  -  t r ( S E“ 1 )  -  log  |S|  -  q  (5) 

where  q  is  the  dimension  in  the  data  (q  =  3  for  3-component  data),  and  tr  = 
trace  of  the  matrix.  The  ML-estimator  is  preferred  mainly  because  it  is 
possible  to  test  the  validity  of  the  models  under  the  Gaussian  assumptions 
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and  white  noise.  Let  F  denote  the  minimum  in  equation  (5).  From  the  general 

properties  of  ML-estimators  it  follows  that  (N-l)F  is  asymptotically 

2 

distributed  as  chi-square  (x  )  with  degrees  of  freedom  equal  to  6  minus  the 
unknown  parameters  and  N  =  (T2  -  +  1). 

2.2  Signal  information  extracted 

The  essential  information  extracted  for  a  given  time  window  is  the 
2 

X  -probabilities  of  P,  SV,  SH,  L  or  R-wave  presence  in  the  3-component 
records  as  a  function  of  time  and  azimuth  as  illustrated  for  P  in  Fig.  1. 

In  plots  like  this,  the  probability  measure  may  be  replaced  by  the  apparent 
velocity  (only  P  and  SV)  because  the  A  matrices  contain  powers  for  the 
dimensions  spanned  by  the  signal.  The  x2-tests  are  mutually  exclusive; 
e.g.,  P  and  SV  cannot  be  triggered  simultaneously.  Finally,  the  estimated 
probabilities,  that  is,  the  max.  value  for  each  time  step,  can  be  used  as  a 
simple  weight  function. 


2.3  Discussion 


Our  approach  to  wavefield  decomposition  or  3-component  record  analysis  is 
in  several  respects  similar  to  those  of  other  investigators,  e.g.,  see 
Christoff ersson  et  al  (1986)  for  references.  A  real  problem  here  is  whether 
further  refinements  of  3-component  record  analysis  are  feasible.  This  would 
entail  use  of  more  complex  particle  motion  models  requiring  more  obser¬ 
vational  degrees  of  freedom,  which  is  equivalent  to  processing  of 
3-component  array  data.  Our  remark  here  is  that  for  the  time  being  we  favor 
extensive  analysis  of  observational  data  rather  than  attempting  further 
methodological  refinements. 

3.  PART  2  -  PRACTICAL  APPLICATIONS 


In  the  subsequent  sections  we  demonstrate  practical  applications  of  our 
3-component  analyzing  technique,  that  is,  wavefield  decomposition  results 
for  event  and  scattering  source  locations,  and  noise  directionality 
studies . 
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3.1  Event  locacions  using  broad-band  records 


The  conventional  approach  to  event  location  is  that  of  using  P  arrival 
times  from  networks  of  stations.  For  arrays  and  single-site  3-component 
stations  bearing  (azimuth)  and  epicentral  distances  estimates  are  needed 
for  event  locations.  Thi  entails  that  relatively  more  information  is 
extracted  from  the  records  as  schematically  illustrated  in  Fig.  1.  The  out¬ 
come  of  an  experiment  of  event  locations  using  the  NORESS  broad-band  recor¬ 
dings  of  6  large  earthquakes  in  the  distance  range  40-155°  are  illustrated 
in  Table  1  (e.g.,  for  details  see  Husebye  et  al,  1986  -  submitted  BSSA) . 


Event  locations  at  local  and  regional  distances 


In  such  cases,  data  analysis  is  tied  to  short  period  recordings,  say 
2-12  Hz,  determining  epicenter  locations  is  with  azimuth  from  estimated 
probabilities  in  the  time/azimuth  plane  and  distance  derived  from  differen¬ 
tial  travel  times  of  phases  like  Pg,  Pn,  Sg,  Sn,  Lg,  etc.  An  instructive 
example  of  3-component  analysis  of  local  events  is  given  in  Fig.  2.  For 
comparison,  semblance  analysis  of  the  corresponding  NOkESS  array  recordings 
per  se  is  also  given.  With  extensive  practical  experience  we  consider  it 
feasible  to  locate  events  using  single-site  3-component  records  within  at 
least  50  km  of  "true"  epicenter  locations  (e.g.,  for  details  see  Ruud  et 
al,  1986,  manuscript  in  preparation). 


3.3  Event  locations  at  teleseismic  distances  -  short  period  records 


In  such  cases  the  distance  estimate  has  to  be  derived  from  apparent  veloci¬ 
ties  alone,  and  thus  large  errors  of  10-15  deg  may  ensue.  However,  prac¬ 
tical  experience  say  in  terms  of  probability/velocity  distributions  of  past 
events  in  the  same  general  area  combined  with  seismicity  knowledge  may 
ensure  reasonable  event  locations  in  such  cases  as  well  -  say  within 
5-10  deg. 
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3»4  Location  of  scattering  sources 


It  is  well  established  that  significant  parts  of  P-  and  S-wave  codas  repre¬ 
sent  scattering  contributions  from  secondary  sources  in  the  lithosphere 
both  at  the  source  and  receiver  ends.  Although  theoretical  treatments  con¬ 
sider  such  contributions  to  be  randomized,  more  prominent  ones  would 
"behave  as  deterministic",  that  is,  being  repeatable.  Since  the  3-component 
analysis  "identifies"  secondary  wavelets,  those  with  Pg-velocities 
(receiver  end)  have  been  located  for  a  suite  of  15  teleseismic  events.  The 
essence  of  these  results  was  that  dominant  scatterers  were  found  in  NE 
quadrant  and  at  distances  of  100-300  km.  An  intriguing  observ  t.on  is  that 
transverse  (SH)  motions  are  occasionally  found  in  the  middle  of  the  P-wave 
train,  entailing  a  depth  to  an  anomalous  body  beneath  NORESS  at  15-20  km, 
but  never  seen  on  all  the  four  3-component  stations  at  A0,  C2,  C4  and  C7 
for  one  event. 


3.5  Noise  directionality 

Practical  3-component  analysis  often  resulted  in  "triggering"  prior  to  P- 
wave  onset.  This  phenomenon  can  be  associated  with  hydroelectric  power  sta¬ 
tion  operation  at  Hunderfossen  (A  ~  73  km,  az  ~  310  deg),  and  Braskereidf oss 
to  the  east  (A  ~  14  km,  az  ~  105  deg).  The  former  is  characterized  by  Pg 
and  Sg  phases  while  the  latter  with  Rg  phases  (vel  ~  2.7  km  s-*)  as 
illustrated  in  Fig.  3.  Noise  directionality  is  occasionally  observed  for 
strong  surf  on  the  coast,  being  of  the  Sg/Lg  type  and  most  prominent  in  the 
1-3  Hz  frequency  band  (e.g.,  see  Ingate  et  al,  BSSA,  1985). 

4.  FINAL  REMARKS  -  CONCLUSION 


The  high-quality  records  from  NORESS  testify  towards  the  source  and  propa¬ 
gation  complexities  of  high-frequency  seismic  waves.  The  key  to  a  proper 
understanding  and  analysis  of  such  phenomena  is  that  of  parameterizing  to 
the  largest  possible  extent  the  information  in  the  event  records.  In  this 
respect,  extracting  polarization  characteristics  of  the  wavefield  is  highly 
significant  as  demonstrated  here. 
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Evenc 

Origin  Time 
h  m  s 

Lac . 
(deg) 

Long. 

(deg) 

H 

(km) 

mb/err 

(deg) 

.  Disc, 
(deg) 

Azi 

(deg) 

Vel. 

(km/ s) 

Region 

6/b-85 

(U 

02.40.12.8 

(02.40.05) 

0.95N 

(0.20N) 

28.43W 

(28.70W) 

10 

(33) 

6.3 

(0.80) 

67.10 

(67.90) 

224.21 

(224.20) 

17.4 

(16.8) 

Cencral  Mid- 
Adancic  Ridge 

29/7-85 

(2) 

07.54.44-3 

(07.54.42) 

36.19N 

(37.54N) 

70.89E 

(71.54E) 

101 

(33) 

6.7 

(1-45) 

42.10 

(43.50) 

93.30 

(94.0) 

13.60 

(11.60) 

Hindu  Kush  Region 

23/8-85 

(3) 

12.41.59.7 

(12.42.02) 

39.42N 
(39 .3  IN) 

75.27E 

(75.39E) 

33 

(33) 

6.4 

(0.14) 

43.9 

(44.0) 

89.0 

(89.0) 

13.30 

(11.50) 

Southern  Xinjang, 

China 

2  1/9-85 
(4) 

01.37.13.8 

(01.37.15) 

17.82N  101. 67W 

(16.94N)(102.61W) 

33 

(33) 

6.3 

(1.26) 

85.21 

(66.40) 

298.59 

(299.0) 

22.4 

(19.30) 

Near  coast  of 

Guerrero,  Mexico 

5/10-85 

(5) 

15.24.02.2 

(15.23.57) 

62.26N  124. 31W 

(61 . 29N) ( 1 24 .94W) 

10 

(33) 

6.5 

(1.01) 

52.49 

(53.50) 

335.88 

(335.10) 

15.0 

(13.5) 

Northwest  Territories, 
Canada 

7/11-85 

(6) 

19.12.29.8 

(19.12.08) 

35.20S  179. 36W 

(35.18SX178.16W) 

33 

(33) 

6.2 

(0.98) 

153.5 

(153.7) 

20.3 

(18.1) 

32.0 

(55.1)) 

East  of  North  Island, 

New  Zealand 

TABLE  1.  Focal  parameters,  taken  from  the  PDE  listings  of  USGS,  for  the  events 
used  in  analysis.  Parameters  in  parentheses  are  those  estimated  from  the  single¬ 
site  3-comp,  broadband  recordings  (NORESS).  The  distance  differences  between  the 
two  sets  of  event  location  estiamtes  are  listed  in  the  mb  column. 


FIG.  1.  3-component  analy¬ 
sis  of  a  Hindu  Kush  earth¬ 
quake  (Event  2  in  Table  1). 
Probabilities  are  plotted 
as  a  function  of  time  (min) 
and  azimuth  (deg).  Window 
length  is  20  sec,  time  and 
azimuth  increments  are  5 
sec  and  1  deg,  respec¬ 
tively.  Original  but 
rotated  broadband  traces 
are  displayed  at  top;  the 
stippled  ones  are  probabi¬ 
lity  filtered.  From  the 
figure  and  results  in  the 
table,  we  have  that  azimuth 
and  distance  (differential 
P,  PP,  S  times)  resolutions 
are  very  good  for  low- 
frequency  body  waves. 
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FIG.  2.  (a)  3-component  probability  plot  for  local  explosion  (A~3  deg)  for 
site  C4;  (b)  semblance  analysis  using  whole  NORESS  with  azimuth  fixed  at 
248  deg,  same  rime  interval.  Note  similarity  in  results,  i.e.,  iden¬ 
tification  of  major  Pn  and  Pg  wavelets. 


FIG.  3.  Noise  directionality  as  illustrated  for  presumed  vibrations  at 
Braskereidf oss  power  station  at  06.47,  27  Nay  1985  (spring  flooding  in 
river  Glomma?).  Semblance  analyzed  used  with  vel  =  2.7  kms-^;  triggering 
occurred  only  at  az~100  deg.  The  associated  beam  wavelets  are  also 
displayed.  Such  effects  are  also  manifested  in  the  NORESS  detection  log, 
and  naturally  generate  Pg-triggering  in  3-component  analysis. 
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OBSERVATIONAL  AND  THEORETICAL  STUDIES 
OF  REGIONAL  WAVE  PROPAGATION 


C.A.  Langston,  A.  Lakhtakia,  V.K.  Varadan,  and  V.V.  Varadan 

Department  of  Geosciences 
Department  of  Engineering  Mechanics 
The  Pennsylvania  State  University 
University  Park,  Pennsylvania  16802 


Introduction 


The  research  program  discussed  here  is  concerned  with  several 
short-period  wave  propagation  topics  which  have  direct  bearing  on 
explosion  discrimination  and  location.  To  repeat  an  obvious  but 
important  point,  determination  of  source  parameters  is  controlled 
by  the  knowledge  one  has  of  the  wave  propagation.  There  is  a 
complete  trade-off  between  source  and  structure  effects  so  one 
has  to  be  known  to  estimate  the  other.  In  our  current  research 
program,  we  have  been  developing  and  using  methods  to  treat  wave 
propagation  at  regional  distances  in  a  variety  of  elastic  earth 
models . 


Analysis  of  Small  Events 

Source  depth  is  a  very  important  parameter  for  the 
discrimination  problem  since  it  is  unlikely,  for  practical 
reasons,  that  explosion  sources  will  ever  be  buried  at  depths 
greater  than  a  kilometer  or  so.  Clear  evidence  of  greater  source 
depth  in  a  seismic  signal  nearly  guarantees  that  the  source  is  an 
earthquake.  We  have  undertaken  an  analysis  of  data  from  a  series 
of  aftershocks  from  the  1968  Meckering  earthquake  of  Western 
Australia  to  try  to  determine  obvious  diagnostic  depth  phases  in 
the  short-period  data.  These  events  are  being  used  for  several 
reasons.  Although  the  main  event  is  interesting  in  its  own 
right,  crustal  structure  between  the  source  area  and  the  WWSSN 
station  MUN  where  the  aftershocks  were  recorded  is  well  known  and 
relatively  homogeneous.  The  area  lies  at  the  edge  of  the  Yilgarn 
block  which  is  part  of  the  Australian  Shield.  Topography  is 
gentle  and  the  crustal  section  is  characterized  by  high  velocity 
crystalline  rocks  with  virtually  no  sedimentary  section.  The 
area  has  been  traversed  by  several  refraction  and  reflection 
surveys  so  that  crustal  wave  velocities  are  well  known  (Mathur, 
1974).  Thus,  this  knowledge  of  the  crustal  structure  can  be 
used  in  an  attempt  to  understand  the  nature  of  the  short-period 
seismogram  at  near-regional  distances  through  the  computation  of 
synthetic  seismograms. 

Figure  1  displays  synthetic  seismograms  of  the  vertical 
component  of  motion  at  distances  of  60  to  100  km  for  a  shallow 
dipping  thrust  fault.  Mundaring  station  was  about  80  to  90  km 
from  the  aftershock  zone.  The  choice  of  mechanism  is  based  on 
the  focal  mechanism  of  the  main  M6 . 8  earthquake  but  is  not 
important  for  the  purposes  of  this  discussion.  The  Meckering 
earthquake  was  caused  by  rupture  along  a  shallow  dipping  thrust 
f ault ( Vogf jord  and  Langston,  1986)  which  ruptured  the  surface  but 
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was  less  than  6  km  deep  at  its  deepest  extent.  Thus,  it  is 
likely  that  aftershocks  were  within  10  km  of  the  surface. 

The  synthetic  seismograms  were  computed  using  a  wavenumber 
integration  technique  (Barker,  1984)  which  yields  the  entire 
seismogram.  The  earth  model,  derived  from  reflection  and 
refraction  profiles,  is  displayed  in  Table  1.  The  synthetic 
seismograms  show  that  the  character  of  the  short-period  waveform 
changes  dramatically  with  source  depth  (Figure  1)  At  2  km 
depth,  the  waveforms  are  simple  showing  what  would  be  interpreted 
as  "P"  and  "S"  wave  arrivals.  The  ”P"  arrival  actually  consists 
of  a  number  of  crustal  reflections,  the  direct  wave,  and  the  Sp 
wave  which  travels  horizontally  under  the  free  surface.  The  "S" 
wave  is  composed  primarily  of  the  fundemental  mode  Rayleigh  wave 
which  is  efficiently  excited  for  this  shallow  depth  source.  This 
phase  should  more  correctly  be  called  the  Rg  phase.  Synthetic 
seismograms  for  a  source  at  10  km  depth  show  that  the  fundemental 
mode  Rayleigh  is  not  well  excited  but  is  included  as  a  minor 
arrivel  in  the  S  wavetrain.  More  importantly,  notice  the 
appearance  of  a  relatively  large  secondary  arrival  after  the 
direct  P.  This  is  the  Sp  wave  which  travels  along  the  free 
surface.  The  increase  in  source  depth  over  the  2  km  case  has 
caused  this  phase  to  move  back  in  time  relative  the  direct  wave. 
Thus,  this  phase  can  be  used  to  directly  estimate  source  depth  by 
measuring  its  arrival  time  relative  to  the  first  arrival. 

Data  from  Meckering  aftershocks  are  available  from  the 
Mundaring  WWSSN  station.  Two  such  events  are  shown  in  Figure  2. 
Notice  that  event  A  shows  a  very  simple  waveform  with  a 
distinctive  Rg  phase  arriving  after  the  S  wave.  This  implies  a 
very  shallow  source  depth  according  the  synthetic  seismogram 
calculation.  The  P  wave  of  event  B  shows  two  distinct  phases 
early  in  the  seismogram  which  suggests  the  P  and  Sp  phases.  This 
event  is  interpreted  to  be  deeper  than  event  A. 

The  study  mentioned  here  is  in  an  early  stage  of  development. 
Other  questions  which  must  be  addressed  are  the  influence  of 
source  mechanism  on  the  relative  amplitudes  of  phases  in  the 
short-period  record,  the  effects  of  perturbations  in  the  velocity 
model,  and  effect  of  anelastic  attenuation.  It  is  probably  not 
possible  to  retrieve  standard  source  parameters  such  as  the  focal 
mechanism  from  such  a  limited  data  set  without  extremely  good 
control  on  crustal  structure.  However,  the  mere  existence  of  an 
easily  observable  secondary  phase  may  be  as  good  a  depth 
estimator  as  pP  is  in  teleseismic  short-period  data.  If  phases 
of  this  type  can  be  extracted  from  the  data,  then  this  technique 
would  offer  an  attractive  and  simple  discriminator  of  explosions 
from  earthquakes. 

In  related  work  we  are  computing  synthetic  seismograms  for 
regional  phases  at  larger  distances.  In  particular,  we  are 
interested  in  the  stucture  effect  on  relative  phase  amplitudes 
which  will  affect  yield  estimates.  Figure  3  shows  vertical 
component  seismograms  recorded  at  Urumqui,  PRC  for  5  large 
Semipalat insk  explosions.  Epicentral  distance  is  950  km  and  the 
seismograms  show  relatively  large  P  phases  compared  to  Pg.  These 
data  were  obtained  from  the  PRC  and  were  recorded  by  the 
broad-band  Soviet  Kirnos  instrument.  Hence  the  spectacular  sweep 
in  frequency  seen  in  the  data  from  1  hz  P  waves  to  30  second 
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Rayleigh  waves.  The  large  P  wave  suggests  that  structure  in  the 
upper-most  mantle  along  this  path  is  similar  to  structure  in  the 
Eastern  United  States  where  velocity  increases  with  depth 
producing  high  amplitude  turning  rays  (Langston,  1982).  The 
relative  excitation  of  these  regional  phases  are  being  explored 
for  Asian  structures. 

Theoretical  Scattering  Problems 

A  part  of  our  work  is  concerned  with  obtaining  accurate 
solutions  to  difficult  wave  scattering  problems.  Using  the 
"T-matrix"  or  extended  boundary  condition  method  of 
Waterman( 1975 )  and  Varadan  and  Varadan ( 1 980 )  we  have  solved  the 
problem  of  SH  and  P-SV  line  sources  located  in  a  waveguide  with  a 
corrogated  free  surface.  The  particular  problem  attacked  in  this 
initial  effort  was  posed  to  examine  Lg  and  Pg  scattering  due  to 
topography  on  and  in  the  crustal  waveguide.  Model  parameters  were 
chosen  to  represent  typical  average  crustal  velocities  and 
topography  seen  in  the  Basin  and  Range  province  of  Western  U.S. 

We  are  currently  examining  these  newly  obtained  solutions  to  look 
at  the  effect  of  topographic  scattering  on  wavetrain  durations 
and  effective  attenuation.  Although  2D  solutions,  these  results 
are  useful  in  examining  the  physics  of  scattering  and  will  be 
incorporated  in  analysis  of  Pg  and  Lg  data  from  explosions  and 
earthquakes  in  the  Western  U.S. 

The  T-matrix  method  is  an  exact  method  insofar  as  there  are 
no  theoretical  approximations  made  in  the  solutions.  Its 
extension  to  simple  3  dimensional  structure  is  straightforward 
since  it  makes  use  of  orthogonality  conditions  in  relating  the 
coefficients  of  the  plane  wave  expansion  of  the  incident  wave 
field  to  the  coefficients  of  the  expansion  of  the  geometry  of  the 
scatterer.  In  this  sense  the  method  is  simpler  than  many  common 
methods  of  computing  plane  layer  responses.  Even  so,  there  is 
probably  a  practical  limitation  to  the  method  when  many 
deterministic  scatterers,  such  as  many  crustal  layers,  are 
considered . 

Figure  5  shows  the  geomtry  of  the  line  source  problem.  An 
elastic  layer  over  halfspace  model  is  assumed  with  a  corrogated 
free  surface.  The  line  source  is  buried  within  the  layer  at  y=0 
and  x=0 .  Figure  5  shows  some  sample  results  of  the  scattering 
calculation.  The  amplitude  curves  are  for  P  and  SV  waves 
radiated  into  the  half space.  Note  that  the  presence  of  the 
corrogated  free  surface  significantly  increases  the  amount  of 
shear  wave  radiated  into  the  halfspace  for  high  angles  of 
incidence.  Surface  displacements  have  also  been  computed  using 
these  solutions  to  examine  the  scattered  field  for  waves  in  the 
crustal  waveguide. 
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Table  1 


Southwestern  Australia  crustal  model  used  in  the  synthetic 
seismogram  calculations (Mathur ,  1974) 


Vp ( km/sec ) 

Va  ( km/sec ) 

Rho ( gm/cc ) 

Qp 

Qa 

Thickness ( km ) 

6. 13 

3.54 

2.78 

1000 

500 

10 

6.70 

3.87 

2.94 

1000 

500 

7 

7.49 

4.32 

3.  10 

1000 

500 

25 

8.39 

4.84 

3.45 

1000 

500 

-- 
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VERT  DEPTH  =  2  km  VERT  DEPTH 
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Hqure  1:  Wavenumber  integration  synthetic  seismograms  for  a  thrust  fault  point  source  at 

*■  m  (ft>  and  10  km  (B!  depth.  The  vertical  component  of  motion  is  shown  as  a  function  of 
range.  Major  phases  are  denoted  on  the  bottom  seismograms. 


Figure  2:  Two  Meckering  aftershocks  recorded  at  the  WWSSN  station 
MUN  (Mundaring,  Australia)  at  a  distance  of  80  to  90  km. 
Aftershock  A  shows  a  characteristic  simple  waveform  with  simple  P 
and  S  arrivals  and  a  very  large  Rg  phase.  Synthetic  seismogram 
calculations  suggest  that,  irregardless  of  source  mechanism, 
sources  at  2  km  depth  and  less  will  display  this  kind  of 
waveform.  Event  B,  however,  shows  a  double  P  arrival  suggesting 
P  and  the  SP  surface  phase.  The  Rg  phase  is  relatively  small, 
compared  to  the  body  wave  phases,  suggesting  a  larger  source 
depth.  The  SP  -  P  time  suggests  a  source  depth  of  5  km. 
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Figure  3:  Vertical  component  of  motion  for  5  Semipalatinsk 
explosions  recorded  at  Urumqi,  PRC.  Distance  is  950  km. 
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Figure  5:  Famzone  P  wave  and  SV  wave  radiated  amplitudes,  in  the 
lower  half-space,  due  to  the  P-wave  line  source  located  at  r  =(0,0). 
These  amplitudes  are  drawn  as  functions  of  6,  the  angle  with  respect 
to  the  -y  axis.  The  surface  Ta  is  specified  by  y=a-bcos(27rx/L)  while 
is  the  plane  y=-d.  The  relevant  parameters  are  vp=6>  vg=3.5  km/sec, 
p=  2.7  gm/cm  for  the  layer,  and  v  =8,  v  =4.5,  p=  3.2  for  the  halfspace. 
Also  a/L  =  d/L  =  1.0  and  o/L  =  0.05.  The  normalized  frequency  is 
(A)  kpL  =  5.0,  fB)  kpL  =  7.5,  (C)  kpL  =  10.0.  The  left  hand  side  shows 
amplitude  for  a  planar  free  surface,  the  right-hand  side  for  free-surface 
corrocation. 
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INTRODUCTION 

A  number  of  recent  studies  have  revealed  that  Lg  measurements  can  provide 
remarkably  stable  and  precise  estimates  of  relative  yield  of  underground  nu¬ 
clear  explosions  (Alexander,  1984;  Nuttli,  1986).  The  two  key  issues  in  yield 
estimation  with  Lg  amplitudes  are  absolute  yield  determination  accuracy  and 
precision  in  relative  yield  determination.  In  order  to  use  Lg  amplitudes  for 
absolute  yield  estimation,  i.e.,  to  use  Lg-magnitude  versus  yield  relations 
for  NTS  explosions  to  estimate  yields  of  explosions  in  another  part  of  the 
world,  it  is  necessary  to  correct  Lg  magnitudes  for  propagation  effects.  Most 
propagation-path  correction  methods,  such  as  that  of  Nuttli  (1986),  assume 
that  only  anelastic  attenuation  and  geometric  spreading  cause  Lg  amplitudes  to 
decay  with  distance.  However,  geologic  blockages  are  also  known  to  profoundly 
affect  Lg  propagation  (Ruzaikin  et  al,  1977)  ,  and  such  effects  are  not  cur¬ 
rently  accounted  for  '.n  Lg  propagation  path  corrections.  The  relative  preci¬ 
sion  of  estimating  yields  with  Lg  amplitudes  has  been  improved  by  making 
measurements  on  digital  seismic  traces.  Digital  prefiltering,  spectral- 
magnitude  measurements,  and  time-domain  measurements  of  window-averages  or 
rms-Lg  amplitudes  has  permitted  more  stable  magnitude  estimates  than  peak-to- 
trough  type  measurements  on  analog  traces. 

The  purpose  of  this  study  has  been  to  investigate  yield-determination 
with  Lg  amplitudes  utilizing  digital  single-channel  and  array  data.  The  three 
primary  objectives  include:  (1)  to  investigate  the  spectral-slope  method  for 
determining  Lg  attenuation,  (2)  to  determine  if  Lg-Q  measurements,  determined 
by  this  method,  are  predictive  of  Lg  amplitude  variations  for  different  paths 
in  Eurasia,  and  (3)  to  develop  amd  study  analysis  techniques  for  using 
digitally  recorded  Lg  waves  for  yield  estimation. 
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Lg  AMPLITUDE  VARIATIONS  IN  EURASIA 


NORSAR  data  for  17  presumed  Soviet  PNEs  and  presumed  weapons  tests  at 
Semipalatinsk  were  used  in  this  study.  Their  locations  are  shown  in  the  top 
map  in  Figure  1.  The  bottom  map  shows  some  of  the  major  geologic  features  in 
western  Eurasia.  Although  western  Eurasia  is  generally  viewed  geologically  as 
a  stable,  laterally  homogeneous  craton,  there  are  a  number  of  major  tectonic 
features,  in  particular  large  sedimentary  basins,  which  may  have  an  important 
effect  on  Lg  propagation. 

Figure  2  shows  a  set  of  incoherent  beam  plots  for  SW-NE  trending  line  of 
events  in  northwestern  Russia.  The  plots  were  made  by  stacking  log-rms  ampli¬ 
tudes  made  on  each  NORSAR  channel  aligned  to  P.  Each  channel  was  prefiltered 
from  0.6  to  3.0  Hz  prior  to  stacking.  Figure  2  shows  that  the  Sn~Lg  wave 
train  is  well  recorded  for  the  first  two  events.  However,  beginning  with 
event  6  (August  11,  1984)  ,  the  Lg  amplitudes  are  sharply  reduced.  As  can  be 
seen  from  Figure  1,  the  Lg  waves  from  events  6,  2A,  14A,  and  13A  must  cross 
the  Timan  ridge,  the  Pechora  sedimentary  basin,  and,  in  the  case  of  event  13A, 
the  northern  Urals.  We  argue  that  these  structures  partially  block  Lg  and 
reduce  its  amplitude  more  than  what  would  be  expected  from  anelastic  attenua¬ 
tion  and  geometric  spreading. 

From  studies  of  incoherent  beams  of  other  events,  we  have  observed  block¬ 
ages  associated  with  the  Pri-Caspian  sedimentary  basin  and  the  Ural  Mountains. 
The  top  map  in  Figure  4  summarizes  our  observations.  We  conclude  that  lateral 
variations  in  the  thickness  of  the  sedimentary  column  above  the  basement  is 
the  dominant  geological  effect  on  Lg  propagation  efficiency  across  western 
Eurasia. 


Lg  SPECTRA  AND  Q— ESTIMATION 

We  have  used  the  spectral-slope  method  described  by  Tang  and  Alexander 

(1985)  to  estimate  frequency-independent  Q  for  Lg.  All  Lg  spectra  were 

computed  for  a  51.2  second  window  starting  at  about  the  3.5  km/sec  group- 

velocity  time.  Spectra  for  each  channel  were  smoothed  and  averaged  over  the 

2 

array.  The  array-average  spectra  were  then  corrected  for  noise  and  an  o> 
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explosion  source  model.  The  high-frequency  decay  of  the  corrected  spectra 
were  then  fit  by  eye  to  the  linear,  frequency-independent  Q  model  of  the  form 

'  109  Ao  -  (Ylf)  £ 

where  A  is  distance  in  km,  U  is  the  group  velocity,  assumed  to  be  3.5  km/sec, 
and  f  is  frequency  in  Hz. 

The  array  average  spectra  for  the  07/09/72  (event  5A)  and  10/04/71  (event 
3A)  PNEs  are  shown  in  Figure  3.  The  bottom  plot  shows  the  two  spectra  on  the 
same  scale  as  the  average  of  the  noise  .pectra  between  the  two  events,  plotted 
as  a  dashed  line.  In  the  top  plot  in  Figure  3,  the  spectrum  of  event  5A  has 
been  increased  at  all  frequencies  by  one  log  amplitude  unit  in  order  to 
separate  it  from  event  3A.  Both  these  events  occurred  in  the  western  Russian 
platform,  their.  Lg  waves  look  very  similar  (see  Figure  2)  and  have  not  under¬ 
gone  any  blockages.  Yet  event  3A  seems  to  have  an  anomalously  low  Q  of  585 
with  y  of  0.17  deg-1,  compared  to  that  of  event  5A.  As  shown  in  top  of  Figure 
3,  this  low  Q  for  event  3A  seems  t*.  have  resulted  from  a  spectral  null  at  2  Hz 
which  resulted  in  the  slope  being  underestimated.  We  have  observed  similar 
nulls  in  spectra  of  other  PNEs  and  these  nulls  appear  to  have  been  produced  by 
source  effects,  such  as  multiple  events.  Thus,  we  conclude  that  event  3A  has 
the  same  Lg  Q  as  event  5A. 

The  map  in  Figure  4  summarizes  our  Q  estimation  results.  We  have  found 
that  Q  ranges  from  775  to  1045  with  corresponding  7s  ranging  from  0.1  to  0.15. 
Local  fluctuations  in  Q  do  not  appear  to  correlate  strongly  with  Lg  amplitude 
fluctuations  or  geology.  Some  of  these  variations  may  be  caused  by  Q  estima¬ 
tion  uncertainties  caused  by  spectral  nulls.  However,  our  preliminary  conclu¬ 
sion  is  that  geologic  blockage  effects  on  Lg  amplitudes  are  not  corrected  by 
Lg  Q  estimates. 


Lg  AMPLITUDE  MEASURB4ENTS  FOR  NTS  EXPLOSIONS 

We  have  also  been  studying  Lg  spectral  magnitudes  for  NTS  events  recorded 
at  the  RSTN  stations.  Our  approach  has  been  to  apply  a  suite  of  narrow 
Gaussian  band  pass  filters  across  the  0  to  5  Hz  band  for  time  windows  on  P 
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noise,  P,  P  coda,  Lg-precursor  P  coda,  Lg,  and  Lg  coda.  We  then  compute  rms 
averages  in  the  windows  in  each  band  and  average  the  log-rms  values  over  a 
selected  frequency  band. 

Figure  5  compares  0.75  to  3.0  Hz  spectral  magnitudes  for  P  and  Lg, 
measured  on  the  vertical  component  at  RSSD,  with  the  scaled  depth  (h/Y1/3)  of 
the  explosion  relative  to  the  scaled  depth  of  the  water  table.  The  15  events 
were  tested  at  one  of  two  yields  below  150  kt.  The  X  symbols  represent  the 
high-yield  events  and  the  O  symbols  represent  the  low-yield  events.  Points  to 
the  right  of  the  vertical  dashed  line  represent  events  detonated  below  the 
water  table. 

The  results  show  the  expected  result  that  events  above  the  water  table 
have  lower  P  and  Lg  amplitudes,  due  to  reduced  coupling,  than  events  of 
comparable  yield  below  the  water  table.  Also,  Figure  5  shows  that  P  and  Lg 
amplitudes  correlate  with  increasing  depth  below  the  water  table,  although  the 
depth  dependence  is  stronger  for  P  than  for  Lg.  In  general,  we  have  found 
that  Lg  magnitudes  have  significantly  less  magnitude  yield  scatter  than  P  wave 
magnitudes  for  events  located  near  or  below  the  water  table. 

The  explosions  we  have  studied  with  RSTN  data  have  yields  in  a  limited 
range  below  150  kt.  In  order  to  extend  the  yield  range,  we  are  collecting  and 
calibrating  additional  digital  data  from  LRSM  recordings  of  the  older,  high- 
yield  U.S.  explosions.  This  data,  which  we  are  putting  into  an  explosion  Lg 
database  at  the  CSS,  will  be  combined  with  the  RSTN  data  assess  the  overall 
yield  estimation  precision  of  digitally-measured  Lg  amplitudes. 
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FIGURE  1  Top;  ISC  or  NEIS  locations  of  the  presented  underground  nuclear 

explosions  used  in  this  study.  Bottom;  Map  showing  the  aajor 
geologic  features  in  western  Eurasia. 
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LOG  RMS  AMPLITUDE 


FIGURE  2 


Lg 


TIME  (SEC) 


Array-stacked  log-rms  amplitude  coda  plots  of  EVEs  located  in 
northwestern  Russia.  Tick  aarks  on  solitude  scale  indicate  0.5 
log-ms  units.  Borisontal  lines  indicate  average  P  background 
noise  levels.  Bach  trace  was  pre— filtered  in  the  0.6  to  3  Hz 
band. 
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LOG  AMPLITUDE 


3.00 


A  =  07/09/72  (SA)  Q  =  876 
B  =  10/04/71  (3A)  Q  =  586 
C  =  AVERAGE  NOISE 


FREQUENCY  (HZ) 

Spectra  of  events  5A  and  3A.  Top  plot  shows  the  spectra  separated 
by  one  log  unit.  Lines  oo  the  top  plot  indicate  the  high- 
frequency  slopes  used  to  estimate  Q.  The  dashed  line  in  both 
plots  indicates  the  average  noise  spectrum  between  the  two  events. 
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FIGOKB  3 


FIGURE  4 
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Nap  suMtar ixing  the  relative  S„  and  Lg  excitation  and  spectral  Q 
■easureaents  for  BOBSAR  recordings  fro«  presided  Soviet  PNEs. 
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STATION  RSSD 

INTEGRATED  9. 75-3.0  HZ  RMS  FOR  P 
SCALED  DEPTH  RELATIVE  WATER  TABLE 


VIGORS  5  Plots  of  spectral  Magnitude  for  P  (top)  and  Lg  (bottoa)  versus  scaled  depth 
with  respect  to  the  water  table  for  15  events  at  VTS.  Points  to  the  rights 
of  the  dashed  line  denote  explosions  below  the  water  table.  Open  circles 
indicate  the  lower  yield  events  and  the  X's  indicate  higher  yield  events. 


THE  TECTONIC  RELEASE  SIGNATURE  AT  REGIONAL  DISTANCES 


Terry  C.  Wallace 
Department  of  Geosciences 
University  of  Arizona 

INTRODUCTION 

One  of  the  main  consequences  of  decreasing  the  yield  limit  in  a  threshold 
test  ban  is  the  increased  reliance  on  regional  distance  seismograms  for 
monitoring  purposes.  Most  schemes  for  determining  the  yield  at  regional 
distances  are  based  on  the  amplitude  of  phases  such  as  Pg  or  Lg.  Here  we 
investigate  the  effect  of  tectonic  release  on  the  character  of  Pn-Pg 
wavetrain,  and  whether  it  introduces  any  bias  in  the  yield  determination. 

There  is  an  apparent  frequency  dependence  of  the  tectonic  release 
signature  at  regional  distances.  At  periods  of  5-15  seconds,  Wallace  et  al. 
(1983)  show  that  the  Pn^  waveforms  from  Pahute  Mesa  explosions  are  sometimes 
strongly  distorted  by  tectonic  release,  which  amounts  to  superimposing  the 
wavefield  from  a  strike-slip  shear  dislocation  on  the  explosion  signature.  On 
the  other  hand,  Alexander  (1980),  Pomeroy  et  al.  (1982)  and  Gupta  and  Blanford 
(1983)  present  results  which  indicate  that  tectonic  release  has  little  effect 
on  3  Hz  Lg  amplitudes.  A  similar  problem  is  rectifying  the  far-field 
representation  of  the  tectonic  release  with  the  observations  of  the  tangential 
component  of  strong-ground  motion.  The  observed  near-field  SH  waves  are  much 
more  complicated  and  significantly  smaller  than  expected. 

We  attempt  to  explain  both  the  frequency  dependence  of  the  tectonic 
release  and  the  near- fie  Id/ tele sei sm ic  imcompatability  with  a  distributed 
shear  dislocation  model.  In  the  near-field,  the  modeled  waveforms  are  highly 
complicated  due  to  interfernece  of  the  various  sources.  At  regional  distances 
the  Pg  arrivals  due  to  the  deviatoric  component  of  the  source  are  smaller  than 
expected  from  a  single  point  source  due  to  the  depth  dependence  of  the  Pg 
excitation. 
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FREQUENCY  DEPENDENCE  OF  TECTONIC  RELEASE 
It  is  a  fairly  straight  forward  process  to  invert  for  the  moment  tensor 
(MT)  of  a  seismic  source  provided  adequate  Green's  functions  can  be 
constructed,  although  for  shallow  explosions  high  frequency  data  is  required 
to  constrain  the  Mz^  components.  The  explosion  moment  is  the  average  of  the 
trace  of  the  moment  tensor,  and  the  deviatoric  component  is  the  remainder.  We 
inverted  the  regional  body  waves  recorded  on  the  broadband  Lawrence  Livermore 
National  Laboratory  network  for  three  Pahute  Mesa  explosions  (CHESHIRE,  FARM 
and  POOL)  in  different  frequency  bands  in  an  attempt  to  isolate  any  frequency 
dependence  of  the  deviatoric  component  of  the  MT. 

At  long  periods  (>2  seconds)  a  two  layered  crust  over  a  halfspace  mantle 
was  a  sufficient  travel  path  structure  to  recover  the  moment  tensor.  Table  1 
shows  the  isotropic  and  major  double  couple  components  of  the  MT.  Also  shown 
is  the  size  of  the  minor  double  couple,  an  indication  of  the  likelihood  that 
the  major  double  couple  actually  represents  tectonic  release.  The  inversion 
was  performed  in  a  linear  fashion,  which  requires  that  all  the  elements  of  the 
MT  have  the  same  time  history.  It  was  assumed  that  the  explosion  had  a 
Helmberger-Hadley  time  history  with  k=15,  B=0.  Although  this  is  a  poor 
representation  for  the  tectonic  release  time  history,  it  probably  is 
sufficient  considering  the  uncertainty  of  the  travel  path  structure. 

TABLE  1:  MOMENT  TENSOR  INVERSION  RESULTS 


Event 

date 

Mex  (dyne-cm) 

Mdc 

Mmdc 

FARM 

.01  x  1024 

.003 

.002 

CHESHIRE 

2-14-76 

.67 

.36 

.07 

POOL 

3-17-76 

.42 

.19 

.12  1 
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Figure  1  shows  the  frequency  dependence  of  the  deviatoric  part  of  the 
moment  tensor.  The  LLNL  records  were  low  pass  filtered  at  different 
frequencies,  as  shown  on  the  abscissa.  The  resulting  moment  tensor  had  the 
long-period  isotropic  component  subtracted  from  it,  and  the  largest  remaining 
eigenvalue  was  assumed  to  represent  the  tectonic  release  (the  direction  of  the 
eigenvectors  was  not  considered).  As  can  be  seen,  the  apparent  size  of  the 
tectonic  release  decreases  with  increasing  frequency. 

It  is  difficult  to  determine  the  casue  of  the  frequency  dependence  of 
tectonic  release.  One  possible  explaination  is  a  strong  depth  dependence  of 
Pg  excitation  by  a  strike-slip  source  (the  orientation  of  the  tectonic 
release).  At  high  frequencies  it  is  possible  that  the  distributed  nature  of 
tectonic  release  results  in  destructive  interference. 

NEAR-FIELD  TECTONIC  RELEASE 

Figure  2  shows  the  three  components  of  ground  motion  from  BOXCAR  recorded 
at  a  distance  of  7.^*  km.  Below  the  observations  are  synthetics  calculated  using 
an  explosion  source  function  and  velocity  structure  determined  by  Barker  et 
al.  (1985).  Wallace  et  al.  (1986)  determined  a  seismic  moment  of  0.86  x  10^ 
dyne-cm  for  the  tectonic  release  from  this  event  on  the  basis  of  teleseismic 
modeling.  The  SH  waves  that  would  be  expected  for  this  size  tectonic  release 
are  shown  below  the  observed  tangential  component.  The  gross  mismatch 
indicates  two  things:  (1)  significant  off  azimuth  energy,  and  (2)  the 

inadequacy  of  a  point  source  representation  for  the  tectonic  release. 

On  the  basis  of  our  teleseismic  modeling  of  the  SH  energy  from  nuclear 
explosions,  we  would  argue  strongly  that  there  is  appreciable  release  of 
tectonic  strain,  and  that  the  tangential  energy  is  not  a  product  of  asymmetry 
of  the  explosion  source  (see  Wallace  et  al.  1985;  1986).  If  the  tectonic 
release  is  purely  a  triggered  fault  motion  as  suggested  by  Aki  and  Tsai 
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(1972),  then  it  would  be  expected  that  the  strong-motion  SH  waves  would  be 
similar  to  those  produced  by  a  similar  sized  earthquake.  Figure  3  compares 
the  SH  waves  from  an  aftershock  of  the  1979  Imperial  Valley  earthquake  ( = 
5.0,  MQ  =  0.7  xlO^  dyne-cm,  recorded  at  10.4  km  distance)  to  the  BOXCAR 
record  and  a  synthetic  for  a  simple  dislocation  model.  Apparently  the  point 
source  representation  is  sufficient  for  the  earthquake,  implying  that  at  least 
some  component  of  the  tectonic  release  is  due  to  a  mechanism  other  than  a 
triggered  fault  motion.  Archambeau  (1972)  proposed  a  mechanism  for  tectonic 
release  which  is  due  to  the  loss  of  strength  of  material  in  a  volume  about  the 
explosion . 

In  an  attempt  to  simulate  the  Archambeau  model  we  constructed  a 
spherically  symmetric,  distributed  source.  Figure  4  illustrates  the 
distributed  model  as  16  point  sources  (of  equal  moment)  about  the  working 
point.  The  elastic  radius,  or  extent  of  shear  stress  drop,  expected  for 
Archambeau' s  model  is  given  by: 

2 

Re  =  Fpa  4,{>}  Ma 

where  F  is  the  ratio  of  the  tectonic  release  moment  to  the  explosion  moment, 
Aa  is  the  stress  drop,  p  and  a  are  the  source  region  density  and  P  wave 
velocity  respectively,  and  { °°}  is  the  static  level  of  the  explosion 
displacement  potential.  For  BOXCAR  (  lOOOkt),  and  the  F  factor  calculated  on 
the  basis  of  the  tectonic  release  moment  of  0.87  x  10^  dyne-cm,  the  elastic 
radius  should  be  on  the  order  of  1.1  km  if  the  stress  drop  is  100  bars.  This 
radius  seems  extrodinarly  large,  and  indeed  would  have  a  significant  effect  on 
a  recording  only  7  km  away.  The  value  of  100  bars  for  stress  drop  is  taken 
from  the  "average"  intraplate  earthquake  (Liu  and  Kanamori,  1983).  In  earlier 
work  (Wallace  et  al.,  1983)  it  was  shown  that  the  sP  waveform  from  BOXCAR 
required  a  much  higher  stress  drop:  300-500  bars.  This  is  consistent  with 
the  work  of  McKeown  and  Dickey  (1969)  who  showed  that  the  aftershock 
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distribution  of  explosions  with  large  F  factors  is  confined  to  a  region  much 
smaller  than  expected  for  an  earthquake  of  comparable  moment.  Using  a  stress 
drop  of  500  bars  the  elastic  radius  for  BOXCAR  is  0.3  km.  We  used  this  radius 
to  define  the  sphere  of  the  distributed  sources. 

Figure  5  shows  a  comparison  of  the  observed  tangential  and  distributed 
source  synthetic  waveforms  for  BOXCAR.  Although  the  match  is  not  perfect,  it 
does  preserve  the  character  of  the  interference  pattern.  Also  shown  in  figure 
5  is  the  long-period  teleseismic  SH  pulse  from  the  distributed  source  '-ompared 
with  that  from  a  point  source.  It  is  apparent  that  the  distributed  source 
model  will  satisfy  both  the  local  data  and  far-field  observations.  The  moment 
that  is  required  in  the  near-field  is  .38  xlO^  dyne-cm,  or  2  times  smaller  than 
the  moment  derived  on  the  basis  of  teleseismic  data. 

CONCLUSIONS  AND  SUGGESTIONS 

It  is  apparent  that  underground  explosions  which  have  a  substantial  component 
of  non-isotropic  radiation  are  releasing  the  tectonic  stress  in  a  volume  about 
the  source.  At  the  high  frequencies  which  are  recorded  at  5-50  km  from  the 
source  the  tectonic  release  can  not  be  modeled  as  a  point  source.  There  is 
little  effect  on  the  estimation  of  explosion  source  functions  from  near-field 
data,  (see  Stump  and  Johnson,  1984)  but  it  makes  the  recovery  of  the 
deviatoric  component  of  the  moment  tensor  difficult.  One  possible  scheme  for 
a  stable  yield  estimation  from  near-field  data  would  be  to  first  use  a 
polarization  filter  and  reject  tangential  motion. 

It  is  much  more  difficult  to  explain  the  frequency  dependence  of  tectonic 
release.  At  long-periods  (>2  seconds)  the  regional  distance  body  waves  can  be 
inverted  for  the  deviatoric  component  of  the  moment  tensor,  and  the  results 
are  consistent  with  the  teleseismic  analysis  of  SH.  At  short  periods,  the 
deviatoric  component  of  Pg  is  smaller  than  expected.  At  least  some  of  this 
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reduction  may  be  due  to  unmodeled  propagation  effects,  but  it  is  still  not 
possible  to  rectify  the  explosion  excitation  of  Pg  with  the  lack  of  double 
couple  excitation.  There  are  at  least  two  possiblitites  which  warrent  further 
study:  (1)  a  frequency  dependent  time  function  for  the  tectonic  release,  and 
(2)  asymme' :y  in  the  explosion  or  tectonic  release. 
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Figure  1:  The  ratio  of  the  deviatoric  to  explosion  moment  for  Pahute  Mesa 
explosions  CHESHTRE,  POOL  and  FARM.  Mex  is  only  determined  at  long-periods. 
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Figure  2:  Comparison  of  strong-ground  motion  from  BOXCAR  observed  at  a 
recording  3ite  7.4  km  away  from  the  source  with  synthetic  seismograms 
calculated  using  the  explosion  model  of  Hartzell  et  al.  (1983).  The 
tangential  model  is  from  the  best  fitting  point  source  double  couple. 
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Figure  3:  Comparison  of  BOXCAR  and  Imperial  Valley  aftershock  with  synthetic 
SH  wavetrain.  The  explosion  and  earthquake  have  approximately  the  same  far- 
field  moment. 


i 


Figure  4:  Schematic  of  the  distributed  source  tectonic  release  modpl 
are  shown  in  solid  dots).  All  sources  are  at  a  distance  of  one 
radius  . 
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Figure  5:  Distributed  source  SH  waves  in  the  near-  and  far-field.  The  strong 
motion  fit  requires  a  total  moment  of  0.38  x  102^  dyne-cm.  The  far-field 
comparison  is  for  a  single  point  source  of  moment  0.38  x  lO2^  amd  16  sources 
of  0.22  x  1023  dyne-cm.  There  is  no  difference  in  waveform  or  amplitude. 
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Lg  MAGNITUDES  OF  SELECTED  SHAGAN  RIVER,  EAST  KAZAKHSTAN  EXPLOSIONS 

Otto  W.  Nuttli 
Saint  Louis  University 


INTRODUCTION 

Amplitudes  of  Lg  waves  of  approximately  1-sec  period  provide  a  stable 
estimate  of  magnitude,  m^(Lg),  and  explosion  yield,  Y,  for  Nevada  Test  Site 
explosions.  When  these  relations  are  applied  to  seven  non-NTS  explosions 
of  announced  yield  in  Colorado,  New  Mexico,  Mississippi,  Algeria  and  East 
Kazakhstan,  the  estimated  yields  differ  from  the  announced  values  by  be¬ 
tween  3  and  35%  (Nuttli,  1986). 

In  the  present  study  m^/Lg)  values  are  determined  for  a  selected  number 
of  Shagan  River,  East  Kazakhstan  explosions.  Figure  1  shows  the  location  of 
the  test  site  and  of  WWSSN  stations  whose  data  are  used.  The  closest  sta¬ 
tion  is  NIL  (Nilore,  Pakistan)  at  approximately  1875  km  and  the  most  distant 
is  KON  (Konberg,  Norway)  at  about  4 375  km. 

DATA  AND  RESULTS 

Maximum  Lg  amplitudes,  as  measured  in  the  group  velocity  window  of  3.2 
to  3.6  km/sec,  are  extrapolated  to  an  epicentral  distance  of  10  km  by  the 
formula  (Nuttli,  1986) 

A  (10  km)  =  A(A  )  (&/10)1/3  sin  ( A /ill .  1  )/sin  (lO/lll.l)  2  exp^f(A-  10 

where  A(10  km)  is  the  hypothetical  lg  amplitude  at  10  km  epicentral  distance 
and  A  (A  )  is  the  observed  amplitude  at  distance  A,  in  kilometers.  The  re¬ 
ference  level  of  the  m^(Lg)  scale  was  set  so  that  m^(Lg)  and  m^(p)  values  are 
numerically  the  same  for  eastern  North  America  earthquakes  (Nuttli,  1973). 
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This  leads  to  the  relation 


mb(Lg)  =  5-0  +  log1Q  [A  (10  km)/ll0] 

where  A  (10  km)  is  the  ground  motion,  in  micrometers,  extrapolated  from  a 
distance  A. 

The  coda-Q,  method  of  Aki  and  Ghouet  (1975)  as  adapted  by  Herrmann 
(1980)  was  used  to  obtain  approximate  values  of  Qq  (Q  at  1-Hz  frequency) 
for  each  source-to- station  path.  Station  corrections  then  were  determined 
for  a  selected  number  of  explosions  by  determining  departures  of  individual 
station  mb(Lg)  values  from  the  average  m^Lg)  values.  These  station  correc¬ 
tions  are  in  the  form  of  correction  to  the  Qq  values.  Adjustments  to  the 
coda-Q  determined  values  of  Qq  varied  between  0  and  6.8 %  for  ten  stations. 

mb(Lg)  values  were  determined  for  62  Shagan  River  explosions  between 
1965  and  1985.  Their  mb(p)  values,  as  determined  by  the  ISC  or  the  USGS, 
varied  between  5.4  and  6.2.  A  linear  least- squares  fit  of  the  mb(Lg)  and 
mb(p)  data  gave 

mb(Lg)  =  0.59?  -  0.392  +  (0.892  -  0.067)  m^P) 

where  the  values  following  the  -  signs  are  standard  deviations.  The  data 
do  not  cover  a  large  range  of  mb(p)  values.  If  a  unit  slope  is  assumed, 

as  would  be  expected  from  the  definition  of  the  mb(Lg)  scale,  the  data 

satify  the  relation 

mb(Lg)  =  (-0.036  -  0.122)  +  mb(p) 

For  water- saturated  NTS  explosions,  the  corresponding  relation  is  (Nuttli, 

1986) 

m^Lg)  =  (0.31  -  0.02)  +  mb(p) 
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Combining  the  last  two  equations  results  in  an  estimated  m^P)  bias  of 
0.35  magnitude  units  between  NTS  and  Shagan  River,  assuming  that  the  same 
yield  explosion  in  water- saturated  rock  at  NTS  and  Shagan  River  will  pro¬ 
duce  the  same  m^(Lg)  values. 

Figure  2  shows  m^(Lg)  plotted  versus  time,  from  1976  through  1985»  for 
Shagan  River  explosions.  There  appears  to  be  an  upper  limit,  of  m,  (Lg) 
about  6.1  to  6.2,  for  explosions  after  autumn  1979.  Between  1976  and 
autumn  1979  there  is  a  suggestion  of  a  somewhat  smaller  upper  limit, 
between  about  5*9  to  6.0. 

Table  1  gives  estimates  of  the  yields  of  the  Shagan  River  explosions 
studied,  assuming  that  the  NTS-derived  relations  between  m^Lg)  and  explosion 
yield  for  water- saturated  rocks  are  applicable.  Two  estimates  based  on 
rn^Lg)  values  are  given.  The  first  was  obtained  by  fitting  the  NTS  data 
by  a  quadratic  equation,  and  the  second  by  a  linear  equation  (Nuttli,  1986). 
Also  included  in  Table  1  are  yield  estimates  by  Dahlman  and  Israelson  (1977) 
and  upper-bound  yield  estimates  by  Sykes  and  Cifuentes  (1984). 

For  seven  of  the  fifteen  explosions  for  which  Dahlman  and  Israelson  (1977) 
gave  estimated  yields,  the  m^(Lg)  estimates  are  almost  the  same  as  theirs. 

For  the  remaining  eight  their  estimated  yield  values  are  larger,  in  two  cases 
by  as  much  as  a  factor  of  three.  In  twenty  of  twenty- one  cases  the  Sykes 
and  Cifuentes  (1984)  upper-bound  yield  estimates  are  smaller  than  those 
obtained  by  m,(Lg).  For  the  seven  explosions  of  largest  estimated  yield  by 
Sykes  and  Cifuentes  (1984),  the  m^(Lg)  values  are  within  40%  of  their  values. 
The  percentage  differences  in  estimated  yields  get  larger  for  the  smaller 
explosions. 

Nineteen  explosions  since  April  1976  have  m^(Lg)  values  exceeding  6.00, 
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corresponding  to  an  estimated  yield  of  150  kt.  None  have  m^(Lg)  values 
exceeding  6.21,  corresponding  to  an  estimated  yield  of  300  kt.  The  January 
15,  1965  cratered  explosion  had  an  m^(Lg)-estimated  yield  of  103  to  109  kt, 
compared  to  an  announced  value  of  125  kt  (Marshall  _et  al,  1979). 


CONCLUSIONS 

Lg-wave  amplitudes  recorded  at  WVSSN  stations  in  southern  Asia  and 
Scandinavia  can  be  used  to  determine  m^(Lg)  values  for  Shagan  River  explosions 
of  m,  (p)  i.  5.5.  The  Scandinavian  stations,  at  larger  epicentral  distances, 
have  Lg  amplitudes  five  to  ten  times  larger  than  the  southern  Asian  stations, 
which  indicate  large  differences  in  Q  for  the  two  paths. 

For  well- determined  m,  (Lg)  values,  taken  to  be  those  based  on  the  data  of 
five  or  more  stations,  the  largest  occurred  for  the  December  2?,  1981  event, 
which  had  an  rn^Lg)  of  6.15*  The  corresponding  yield  estimate  is  about  25O  kt. 

The  rn^(p)  magnitude  bias  between  Shagan  River  and  Nevada  Test  Sites,  as 
determined  by  m^(Lg)  values,  is  0.35  units.  This  assumes  that  explosions  of 
the  same  yield  at  the  two  sites  detonated  in  similar  material  would  produce 
identical  m^(Lg)  values.  This  assumption,  when  made  for  explosions  of  an¬ 
nounced  yields  at  five  widely  disparate  test  sites,  gives  estimated  yields 
within  3  to  35%  of  the  announced  values. 
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TABLE  1 


YTELD  ESTIMATES  OF  EAST  KAZAKHSTAN  EXPLOSIONS 

_  Estimated  explosion  yield,  in  kilotons 


Date 

m,(Lg) 

quadratic  eq. 

linear  eq. 

Dahlman  and 

Sykes  and  Ci- 

this  study 

this  study 

Israelson  (1977) 

fue.ntes  (1984)* 

01/15/65** 

5. 87 

103 

109 

110 

— 

02/10/72 

5-55 

42 

42 

43 

— 

11/02/72 

6.04 

171 

183 

350 

— 

12/10/72 

6.09 

200 

212 

620 

— 

07/23/73 

6.13 

225 

239 

420 

— 

12/14/73 

5.87 

103 

109 

150 

— 

05/31/74 

5-68 

60 

62 

140 

— 

10/16/74 

5-26 

20 

17 

43 

— 

12/27/74 

5-69 

62 

64 

51 

— 

04/27/75 

5.47 

34 

33 

60 

— 

10/29/75 

5.45 

33 

31 

90 

— 

12/25/75 

5.83 

92 

97 

90 

— 

04/21/76 

5-19 

17 

14 

20 

— 

06/09/76 

5-27 

20 

18 

25 

— 

07/04/76 

5-90 

113 

120 

90 

— 

08/28/76 

5.60 

48 

49 

— 

— 

11/23/76 

5. 86 

100 

106 

— 

— 

12/07/76 

5-71 

65 

68 

— 

— 

05/29/77 

5.58 

45 

46 

— 

— 

06/29/77 

5-15 

15 

13 

— 

— 

09/05/77 

5.51 

38 

37 

— 

— 

11/30/77 

5-71 

65 

68 

— 

06/11/78 

5-75 

74 

76 

— 

— 

07/05/7 8 

5-67 

59 

60 

— 

— 

08/29/78 

5.80 

85 

89 

— 

61 

09/15/78 

5-87 

104 

109 

— 

59 

11/04/78 

5.57 

45 

44 

— 

10 

11/29/78 

6.01 

156 

167 

— 

69 

06/23/79 

5-92 

120 

127 

— 

186 

07/07/79 

5-87 

103 

109 

08/04/79 

6.01 

156 

167 

— 

146 

08/1 8/79 

6.03 

166 

177 

— 

152 

10/28/79 

6 . 06 

182 

194 

— 

67 
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TABLE  1  (Concluded) 

—  Estimated  explosion  yield,  in  kilotons  - 


Date 

mb(Lg) 

quadratic  eq. 
this  study 

linear  eq . 
this  study 

Dahlman  and 
Israelson  (1977) 

Sykes  and  Ci¬ 
fuentes  (1984)* 

12/02/79 

6.05 

176 

188 

— 

70 

12/23/79 

6.12 

219 

232 

— 

152 

06/12/80 

5.74 

71 

74 

— 

06/29/80 

5-71 

65 

68 

— 

19 

09/14/80 

6.09 

200 

212 

— 

184 

10/12/80 

5-92 

120 

127 

— 

48 

12/14/80 

5.92 

120 

127 

— 

5? 

12/27/80 

6.00 

152 

162 

— 

39 

03/29/81 

5.45 

32 

31 

— 

04/22/81 

5.97 

139 

148 

— 

57 

09/13/81 

6.10 

206 

219 

— 

94 

10/18/81 

6.09 

200 

212 

— 

82 

12/27/81 

6.15 

241 

254 

— 

210 

04/25/82 

6.13 

226 

239 

— 

105 

07/04/82 

6.14 

234 

24  7 

— 

192 

12/05/82 

6.21 

293 

305 

— 

_  __ 

10/06/83 

5.93 

126 

131 

— 

—  —  — 

IO/26/83 

6.10 

206 

219 

— 

02/19/84 

5.74 

71 

74 

— 

—  —  _ 

03/07/84 

5.68 

60 

62 

— 

_  _  w 

03/29/84 

5. 97 

139 

148 

_ . 

04/25/84 

5.86 

100 

106 

— 

_ 

05/ 26/84 

6.07 

188 

200 

_ _ 

___ 

07/14/84 

6.17 

257 

270 

— 

— —  ^ 

10/27/84 

6.10 

206 

219 

— _ 

12/02/84 

5-97 

139 

148 

— _ 

_  _  __ 

12/16/84 

6.08 

194 

206 

_ 

___ 

12/28/84 

6,13 

226 

239 

_ 

_ _ 

02/10/85 

5-98 

144 

152 

— 

_ 

*  Sykes  and  Cifuentes  (1984)  considered  these  upper-bound  estimates. 

**  Marshall  et  al  (1979)  state  that  the  announced  yield  was  125  kt,  in  a 
cratered,  water- saturated  sandstone. 


164 


KEV 


165 


Figure  1.  Location  of  the  Shagan  River  Test  Site  and  of  the  WWSSN  stations  whose  data  were  used. 
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Figure  2.  mb(Lg)  values  of  Shagan  River  Test  Site  explosions  from  1976  to  1985. 


ANALYSIS  OF  TELESEISMIC  P  WAVE  AMPLITUDE  AND  CODA 


VARIATIONS  FOR  UNDERGROUND  EXPLOSIONS 
THORNE  LAY,  CHRISTOPHER  LYNNES,  JOAN  WELC 
DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
UNIVERSITY  OF  MICHIGAN 

A  large  data  set  of  more  than  2200  teleseismic  P  wave  recordings  from  59 
NTS  events  is  being  analyzed  with  the  intent  of  establishing  the  near-source 
contribution  to  observed  amplitude,  travel  time  and  coda  variations.  The  data 
have  been  processed  uniformly,  enabling  one-to-one  correlations  between  magni¬ 
tude  and  travel  time  anomalies,  as  in  Figure  1.  Pahute  Mesa  observations  have 
a  stronger  correlation  coefficient  (0.302  for  individual  observations,  0.364  for 
station  averages)  than  for  Yucca  Flat  events  (0.191  for  individual  observations, 
0.207  for  station  averages),  supporting  a  focussing/defocussing  interpretation 
of  the  Pahute  Mesa  amplitude  pattern.  However,  much  of  the  travel  time  pattern 
originates  far  from  the  source  region,  given  a  strong  correlation  (ccc=0.946) 
found  between  station-averaged  residuals  for  the  two  subsites.  The  station- 
averaged  amplitude  and  travel  time  patterns  for  the  two  sites  are  shown  in 
stereographic  projections  in  Figure  2.  The  differenced  patterns  (Pahute  Mesa  - 
Yucca  Flat)  are  shown  on  the  right,  indicating  the  near-source  component. 

Note  that  the  differenced  patterns  retain  slowly  varying  components,  which  we 
are  attempting  to  explain  deterministically  using  backprojection  and  thin  lens 
modeling  procedures.  Figure  3  shows  an  example  of  backpro jecting  the  anomalies 
to  surfaces  at  different  depths  beneath  NTS.  This  procedure  is  used  to  constrain 
the  spatial  dimensions  of  the  upper  mantle  heterogeneity. 

The  complete  waveform  information  in  the  first  15  sec  of  the  P  arrivals  for 
the  entire  data  set  is  being  analyzed  as  well,  using  statistical  and  deterministic 
procedures.  The  statistical  procedure  involves  computing  the  energy  temporal 
centroid  of  the  P  wave  envelope  for  each  station,  applying  corrections  to  the 
centroid  accounting  for  station  and  source  effects,  and  averaging  the  resulting 
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measurements  for  each  event.  These  average  complexity  (energy  flux). anomalies 
are  then  correlated  with  known  event  parameters  as  in  Figure  4,  in  order  to 
empirically  establish  the  mechanisms  affecting  the  waveforms.  Correlations  with 
position  in  the  test  site  for  Pahute  Mesa  events,  and  with  magnitude  (and  hence 
pp  delay  and  burial  depth)  for  Novaya  Zemlya  events  have  been  detected.  Inspec¬ 
tion  of  the  azimuthal  patterns  of  the  individual  complexity  measurements  shows 
slowly  varying  patterns  for  Pahute  Mesa  that  will  further  constrain  the  causes 
of  the  waveform  variations.  Analysis  of  the  frequency  dependence  of  the  coda 
variations  has  shown  that  some  events  in  Pahute  Mesa  have  azimuthal  patterns  of 
late  high  frequency  arrivals,  possibly  produced  by  early  aftershocks.  Figure  5 
shows  observations  for  GREELEY,  for  which  stations  in  the  loop  directions  of 
the  tectonic  release  radiation  pattern  have  late  high  frequency  energy.  Figure  6 
shows  examples  of  the  azimuthal  patterns  of  the  differences  between  the  envelope 
centroid  times  for  high  and  low  frequency  passbands  for  several  NTS  events.  The 
slowly  varying  patterns,  together  with  the  differences  in  pattern  from  event  to 
event  demonstrate  that  near  source  effects  are  responsible.  We  are  attempting 
to  establish  whether  these  variations  are  due  to  tectonic  release  or  whether 
they  are  related  to  the  upper  mantle  heterogeneity  that  produces  the  first 
arrival  amplitude  patterns. 

An  attempt  to  provide  a  deterministic  interpretation  of  some  of  the  P  wave 
coda  in  the  first  15  sec  of  the  teleseismic  waveforms  is  being  made  using  slant 
stacking  procedures.  The  basic  idea  is  to  seek  scatterer  locations  in  the 
near-source  vicinity.  Energy  ratio  coherency  measures  and  semblance  are  being 
used  to  appraise  the  slant  stacking  output.  Station  and  event  weighting  filters 
have  been  developed  to  reduce  the  influence  of  correlated  noise  in  the  data, 
and  the  full  1300  waveform  data  set  for  Pahute  Mesa  events  is  being  utilized. 
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Figure  2a. 


Figure  2b. 


igure  2.  a)  Lower  hemisphere,  equal  area  projections  showing  the  average 

station  magnitude  and  travel  time  anomalies  for  Pahute  Mesa  events  (left) 
and  Yucca  Flat  events  (right).  Note  the  shift  of  the  amplitude  patterns 
between  the  two  sites  and  the  mutual  shift  relative  to  the  travel  time 
patterns.  The  outer  perimeter  is  for  a  45°  take-off  angle  for  a  reference 
velocity  of  7.8  km/sec.  b)  Similar  projections  of  the  difference  in 
station  anomalies  obtained  by  subtracting  the  Yucca  Flat  values  from  those 
for  Pahute  Mesa. 
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figure  3.  ^ac'-nro j  ections  onto  surfaces  below  NTS  for  amplitude  (left)  and 

travel  time  (right)  anomalies  for  both  Pahute  Mesa  and  Yucca  Plat  events 
The  surfaces  have  been  bin-averaged  to  smooth  the  pattern  for  the  2000 
observations.  Note  that  a  depth  of  greater  than  35  km  is  suggested  for 
the  origin  of  the  common  amplitude  and  travel  time  features  for  the  two 
test  si tes . 
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figure  4a. 
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Figure  4b. 

Figure  4.  a)  Regression  of  average  event  complexity  factors  for  Pahute  Mesa 
events  on  ,  event  magnitude  anomaly,  distance  from  the  center  of  the 
mesa,  and  relative  tectonic  release  F-factor.  b)  Regression  of  average 
event  complexity  factors  for  Northern  Novaya  Zemlya  events  on  m  and 
F-factor. 
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pigure  5.  Example  of  frequency  dependence  of  the  coda  envelopes  for 

observations  for  GREELEY.  The  top  row  is  the  original  trace  for  three 
WV/SSN  stations.  The  middle  row  is  the  narrowband  filtered  trace  with 
passband  O.U  -  0.8Hz,  and  the  bottom  row  is  a  higher  passband  of 
1.1  -  1.5Hz.  The  instrument  response  has  been  deconvolved  in  all  the 
filtered  traces.  Note  the  late,  high  frequency  pulse  for  ESK  and  BOG, 
which  are  in  loop  directions  of  the  tectonic  release  P  radiation  pattern 
ATL  is  along  the  tectonic  release  P  radiation  node.  173 


COMPARISON  OF  MATCHED-FILTERS  TO  TIME-VARiAB LE  FILTERS  FOR  AMPLITUDE 
DETERMINATION  OF  NORMAL  MODES  IN  SURFACE  WAVE  ANALYSTS 


DAVID  R.  RUSSELL 
HORNG-JYE  HWANG 

ST.  LOUIS  UNIVERSITY 


A  currently  accepted  method  for  amplitude  determination  of  normal  modes 
is  phase-matched  filtering  (Stevens,  1986).  The  phase-matched  filter  (PMF) 
involves  centering  a  zero-lag  symmetric  window  in  the  time  domain  on  the 
"pseudo-autocorrelation"  function,  found  by  removing  the  phase  of  a  particular 
mode  of  interest  (Herrin  and  Goforth,  1977).  It  can  be  shown  (Jenkins  and  Watts, 
1968),  that  the  bias  in  the  amplitude  spectrum  caused  from  time-windowing  an 
autocorrelation  function  is  proportional  to  the  second  derivative  of  the  ampli¬ 
tude  spectrum: 

=  (1) 

where 

B  =  bias 

M  =  one-sided  width  of  symmetric  time-domain  window 
A  "  =  second  derivative  of  amplitude  spectrum 
w  =  angular  frequency 

In  this  case  it  is  assumed  that  a  cosine  window  is  used;  other  windows  will  change 
the  coefficient  multiplying  the  second  derivative.  This  bias  can  be  significant  in 
the  vicinity  of  band-edges  of  the  spectrum  or  in  any  region  of  high  curvature.  It 
can  be  reduced  by  extending  the  window  width  M ,  but  this  may  introduce 
extraneous  noise  (or  signals)  into  the  spectrum. 

To  reduce  the  problem  of  bias  and  noise,  the  following  time-variable  filter  is 
constructed: 


f(t)=±/F(U)[S(t  +^L )-S(t  )]c ost-^V-'ctsj 


CJ  I 


4  a 


(2) 


where 
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F(u)  =  pseudo-autocorrelation  function 
S(  )  =  unit  step  function 
a  =  M/T 

M  =  one-sided  width  of  cosine  window  at  period  T 
Taking  the  Fourier  transform  of  /(f)  gives  the  amplitude  spectrum  of  the  mode 
of  interest.  The  above  integral  consists  of  a  sum  of  symmetrically  windowed 
Fourier  harmonics  about  zero-lag  in  the  time  domain.  In  the  traditional  sense  of 
a  TVF,  the  windows  are  symmetric  about  the  group  velocity  peaks  found  from 
multiple  filter  analysis  (Dziewonski  and  Hales,  1972,  p.  74).  Equation  (2)  is  simi¬ 
lar  to  the  traditional  TVF  after  residual  dispersion  is  applied  (Dziewonski  and 
Hales,  1972,  p.  76),  except  that  the  basis  for  this  method  is  phase-matched 
filtering  instead  of  multiple  filter  analysis. 

To  test  the  method,  a  normal  mode  synthetic  seismogram  was  constructed 
at  a  distance  of  6000  km  from  the  source.  The  signal  was  lagged  behind  itself  by 
300  seconds  (to  simulate  multiple  events)  and  superimposed  on  the  original  sig¬ 
nal.  Gaussian  noise  was  then  added  to  the  superimposed  signal  to  give  the  effect 
of  random  noise  (figure  1).  The  noisy  seismogram  was  then  phase-matched 
filtered  about  the  original  signal  and  both  PMF  and  TVF  analysis  was  performed. 
Figure  2  shows  the  effect  of  using  a  M  =  300  second  window  for  the  PMF  and  an 
M =3T  family  of  windows  for  the  TVF.  Due  to  the  wide  window  width,  neither  of 
the  spectra  show  significant  bias,  but  the  TVF  has  a  superior  signal  to  noise  level. 
Figure  3  demonstrates  the  effect  of  narrowing  the  window.  For  the  PMF,  A/=75 
seconds.  The  TVF  still  has  jW  =37',  but  the  maximum  window  width  is  not  allowed 
to  exceed  73  seconds.  At  long  periods  (>  60  seconds),  both  methods  exhibit  bias 
in  the  amplitude  spectra.  However,  the  TVF  is  still  superior  in  reducing  short 
period  noise. 

Another  test  was  made  with  an  earthquake  source  mechanism  at  a  distance 
of  1000  km.  The  fundamental  mode  and  a  single  higher  mode  was  constructed 
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AMP  CCM-SEC)  AMP  CCM-SECD 


Figure  5.  Heavy  lines:  theoretical  fundamental  mode  of  the  earthquake. 


Top:  PMF  with  M=75  seconds. 
Bottom:  TVF  with  M=3T. 
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to  see  the  effect  of  the  two  techniques  on  isolating  the  fundamental  mode  (figure 
4).  In  figure  5,  the  PMF  had  to  have  a  maximum  window  width  of  M  =  75  seconds 
to  isolate  the  fundamental  mode.  Observe  the  bias  at  periods  greater  than  60 
seconds  and  at  20  seconds,  illustrating  the  effect  of  curvature  on  the  windowed 
amplitude  spectrum  (equation  1).  The  TVF  was  constructed  with  M  =  3T  and  no 
long  period  constraints.  This  was  possible  since  the  higher  mode  was  limited  to 
the  shorter  periods.  The  TVF  amplitude  spectrum  shows  little  difference  from 
the  theoretical  fundamental  mode. 

Other  tests  on  synthetics  show  that  it  is  possible  to  reduce  the  width  of  the 
TVF  window  to  M=2T  without  significant  disortion.  However,  these  tests  were 
made  on  given  seismograms,  and  more  analysis  is  needed  to  generalize  the  win¬ 
dow  width  as  a  function  of  curvature  of  amplitude  spectra. 
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0  610  1220  1830  2HH0  3050 


PERIOD  CSECD 

Figure  1.  Synthetic  normal  mode  explosion  recorded  at  6000  km. 

Top:  signal  lagged  and  superimposed  at  300  seconds. 
Center:  Gaussian  noise  added  to  signal. 

Bottom:  Amplitude  spectrum  of  noisy  seismogram. 
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Figure  2.  Heavy  lines:  theoretical  fundamental  mode  of  the  explosion. 
Top:  PMF  with  M=300  seconds. 

Bottom:  TVF  with  M=3T. 
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Figure  3.  Heavy  lines:  theoretical  fundamental  mode  of  the  explosion 
Top:  PMF  with  M=75  seconds. 

Bottom:  TVF  with  M=3T;  window  limited  to  75  seconds  and  less. 
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ANELASTIC  PROPERTIES  OF  THE  CRUST  AND  UPPER  MANTLE 
IN  STABLE  AND  TECTONICALLY  ACTIVE  REGIONS 

H.J.  Hwang  and  B.J.  Mitchell 
Department  of  Earth  and  Atmospheric  Sciences 
Saint  Louis  University 

INTRODUCTION 

An  understanding  of  the  attenuation  rate  of  seismic  waves  with  dis¬ 
tance  is  fundamental  to  our  capability  to  determine  magnitudes  of  earth¬ 
quakes  and  yields  of  nuclear  events.  It  is  also  pertinent  to  detection 
studies,  because  seismic  wave  attenuation  rates  will  govern  the  maximum 
distance  to  which  seismic  waves  may  be  reliably  recorded.  We  have 
determined  the  attenuation  rate  of  seismic  surface  waves  in  several  dif¬ 
ferent  regions  of  the  Earth,  some  of  which  are  relatively  stable,  and 
some  of  which  are  tectonically  active,  in  order  to  determine  the  manner 
in  which  Q  values  and  the  frequency-dependence  of  those  values,  in  the 
crust  and  upper  mantle  vary  regionally.  At  shorter  periods,  our  results 
pertain  to  properties  of  the  crust,  and  are  therefore  important  to  an 
understanding  of  Lg  propagation.  At  longer  periods,  they  provide  infor¬ 
mation  on  the  upper  mantle;  they  therefore  should  contribute  to  our 
ability  to  determine  yields  using  teleseismic  body  waves. 

DATA  ANALYSIS 

Most  determinations  of  seismic  wave  attenuation  are  characterized 
by  large  uncertainties  which  limit  our  ability  to  study  details  of  the 
distribution  of  the  anela3tic  properties,  or  Q,  of  the  Earth.  Hwang  and 
Mitchell  (1986)  recently  showed  that  In  two-station  methods,  measured 
attenuation  coefficients,  at  long  periods,  are  extremely  sensitive  to 
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noise  contamination,  and  they  developed  a  method  which  combines  modal 
isolation  and  frequency-domain  Wiener  filtering,  which  appears  to  pro¬ 
vide  more  stable  and  reliable  determinations  than  methods  previously 
developed. 

New  attenuation  data,  as  well  as  phase  and  group  velocity  informa¬ 
tion,  have  been  obtained  for  eastern  South  America,  western  South  Amer¬ 
ica,  the  Indian  shield,  and  the  Himalayas  (Figure  1).  These  include 
both  two-station  data,  as  well  as  single-station  data,  which  have  been 
combined  to  yield  dispersion  and  attenuation  coefficient  curves  which 
span  the  range  between  a  few  seconds  and  50-100  seconds,  depending  upon 
region.  An  example  of  the  attenuation  coefficient  data  for  fundamental- 
and  first  higher-mode  Rayleigh  waves  in  the  stable  portion  of  eastern 
South  America  appears  in  Figure  2. 

Our  previous  determinations  of  attenuation  coefficients  in  stable 
regions  were  restricted  to  shorter  periods,  30-50s  at  the  long-period 
end,  and  were  accompanied  by  larger  standard  deviations.  As  a  conse¬ 
quence,  our  models  were  very  poorly  constrained  at  larger  depths 
corresponding  to  the  lower  crust  and  upper  mantle. 

VELOCITY  INVERSIONS 

The  phase  and  group  velocity  data  were  inverted  using  a  differen¬ 
tial  inversion  program,  written  by  D.R.  Russell,  which  minimizes  the 
magnitude  of  the  error  vector  of  the  solution  plus  the  differences 
between  adjacent  solution  elements  (Twomey.  1977).  Ir.  addition  to  the 
new  data  described  above,  older  dispersion  data  for  eastern  and  western 
North  America  were  also  reinverted  using  the  new  inversion  scheme.  Par- 
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tial  derivatives  of  phase  velocity  with  respect  to  layer  parameters  for 
these  models  were  used  in  the  Q  inversion  discussed  in  the  following 
section.  The  dispersion  data  for  all  of  the  stable  regions  could  be 
explained  by  models  which  do  not  include  an  upper  mantle  low-velocity 
zone.  Because  of  the  limited  resolving  power  of  the  data  at  upper  man¬ 
tle  depths,  however,  a  minor  low-velocity  zone  cannot  be  excluded  on  the 
basis  of  the  surface  wave  data.  Low- velocity  zones  are  present,  with 
varying  degrees  of  severity  and  resolvability  in  all  of  the  tectonically 
active  regions.  In  regions  where  both  Rayleigh  wave  and  Love  wave  data 
are  available,  there  is  no  need  to  invoke  anisotropy  of  elastic  proper¬ 
ties  to  explain  the  data,  at  least  to  depths  of  about  150  km. 

Q.  INVERSIONS 

The  attenuation  coefficient  data  were  inverted  using  methods  which 
have  been  discussed  in  several  previous  studies  (e.g.  Mitchell,  1980), 
using  the  formulation  of  Anderson  et.  ji.  (1965).  A  model  resulting 
from  the  inversion  of  data  in  eastern  South  America  is  shown  in  Figure 
2.  In  this  stable  region,  as  well  as  others,  we  find  shear  wave  Q 
values  in  the  lower  crust  and  upper  mantle  which  are  lower  than  those 
previously  found  in  stable  regions.  We  attribute  this  difference  to 
improved  data  quality  at  longer  periods,  compared  to  that  of  our  earlier 
studies. 

As  a  rule,  shear  wave  Q  values  in  tectonically  active  regions  are 
significantly  lower  than  those  in  stable  regions,  throughout  the  crust 
and  upper  mantle.  There  is,  however,  a  large  variability  in  the  Q 
models,  even  among  those  for  stable  regions.  In  particular,  shear  wave 
Q  values  at  crustal  depths  beneath  the  Indian  shield  are  significantly 
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lower  that  those  of  other  shields.  These  lower  values  may  be  related  to 
higher  than  normal  reduced  heat  flow  values  which  have  been  determined 
in  the  Indian  shield  (Rao  et.  ai.,  1976).  A  comparison  of  shear  wave  Q 
models  for  several  regions  appears  in  Figure  3« 

FREQUENCY  DEPENDENCE  QE  Q 

It  is  possible  to  study  the  frequency  dependence  of  shear  wave  Q  in 
the  crust,  over  a  limited  period  range,  by  comparing  the  theoretical  Q 
values  of  1-Hz  Lg  waves  produced  by  the  velocity  and  Q  models  described 
above,  with  observed  Q  values  of  Lg.  Theoretical  values  of  Lg  Q  can  be 
obtained  by  computing  short-period  synthetic  seismograms  at  several  dis¬ 
tances,  measuring  amplitudes  of  Lg  at  each  distance,  and  computing  the  Q 
values  needed  to  explain  the  amplitude  fall-off  with  distance.  In  this 
way  we  obtain  a  Q  value  for  Lg  in  much  the  same  way  that  observers 
obtain  it  from  real  seismograms. 

This  process  results  in  theoretical  Q  values  which  are  very  similar 
to  observed  values  in  tectonically  active  regions.  Our  Q  models  for 
those  regions  can  therefore  be  independent  of  frequency  and  satisfac¬ 
torily  explain  observed  surface  wave  attenuation  over  a  broad  period 
range  ( 1  s  to  several  10's  of  seconds).  In  stable  regions,  however,  the 
models  we  obtain  predict  Lg  Q  values  which  are  lower  than  those 
observed.  This  result  indicates  that  shear  wave  Q  is  frequency- 
dependent  in  those  regions,  being  higher  at  a  period  of  1  s  than  at 
longer  periods. 
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Figure  2.  Shear  wave  Q  model  for  eastern  South  America  (upper  left), 

along  with  resolving  kernels  (upper  right).  Comparison  of  observed 
and  theoretical  attenuation  coefficient  values  for  fundamental-mode 
Rayleigh  waves  (bottom  left)  and  first  higher-mode  Rayleigh  waves 
(bottom  right)  are  shown  for  the  model  above. 
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P  WAVE  ATTENUATION,  mb  BIAS  and 
the  THRESHOLD  TEST  BAN  TREATY 

Thomas  C.  Bache,  Steven  R,  Bratt  and  Lolitia  B.  Bache 

Science  Applications  International  Corporation 
10210  Campus  Point  Drive,  San  Diego,  California  92121 

INTRODUCTION 

Concern  about  Soviet  compliance  with  the  150Kt  TTBT  limit  is  motivated  by  the 
large  m^  of  Shagan  River  tests  and  by  the  pattern  of  testing  at  that  site  (Figure 
1)  .  Before  1976  a  broad  range  of  yields  were  tested  at  Shagan  River,  including 
four  with  >  6.0  The  yields  of  Novaya  Zemlya  tests  were  sharply  reduced  after 
31  March  1976.  For  the  first  two  years  of  the  threshold  the  Shagan  River  test 
yields  were  a  factor  of  1.5  -  2.0  less  than  for  previous  tests  at  that  site. 
During  1978  the  yield  increased,  and  many  tests  since  then  have  yields  that  are 
about  a  factor  of  two  larger  than  the  yield  limit  apparently  observed  in  1976- 
1973.  There  are  two  main  competing  hypotheses  to  explain  this  pattern: 

1.  The  bias  between  NTS  tuff/rhyolite  and  Shagan  River  is  large  (0.4- 
0.5)  and  all  Shagan  River  yields  are  150Kt  or  less.  For  some  reason 
the  Soviets  chose  to  restrict  their  testing  to  75-100Kt  during  the 
first  two  years  after  the  limit  went  into  effect. 

2.  The  m^  bias  between  NTS  tuff/rhyolite  and  Shagan  River  is  moderate 

(about  0.2  or  so).  The  events  in  1976-1978  (maximum  m^  -  5.92)  have 

yields  up  to  150Kt.  Events  after  1979  have  yields  up  to  twice  that 
level . 


4> 

E 


7.0  - 

6.8  - 
•.«  - 
•  4  - 

6.2  - 
6.0  - 

5.8  - 

5.6  - 

6.4  - 
5.2  - 
6.0- 

4.6  - 
4.6  - 

4.4  - 


_• _ til- 

•  *  *  X 


%  • 


,5->- 


.  M 
- N 


/. 


•  Shagan  Rtvar 
■  Shagan  fttvar  |N£fS) 
o  Novaya  Zataiya 


1066  1066  1067  1066  1060  1670  1071  1072  1073  1974  1876  1076  1977  1878  1979  1980  1981  1962  1963 


Figure  1  The  m^  are  plotted  versus  time  for  all  large  Novaya  Zemlya  and 
Shagan  River  explosions.  Except  for  the  NEIS  m^,  the  Shagan  River  m^  are  from 
Marshall  et.  al.  (1985)  and  the  Novaya  Zemlya  m,  were  computed  by  Marshall 
(personal  communication)  with  the  same  method.  “The  date  on  which  the  TTBT 
limit  began  to  be  observed  is  marked.  The  horizontal  lines  are  drawn  at  the 
Indicated  m,  . 
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The  major  causes  for  m^  bias  are  differences  in  source  coupling  and  average  path 
attenuation.  Large  source  coupling  differences  exist  (e.g.,  between  saturated 
and  unsaturated  tuffs),  but  given  the  paucity  of  data  from  the  Shagan  River  site, 
we  can  only  make  some  informed  guesses  about  the  contribution  of  source  coupling 
differences  to  bias.  The  average  path  attenuation  differences  between  NTS  and 
E.  Kazakh  seem  easier  to  estimate.  NTS  is  in  a  tectonic  region  known  to  have 
high  attenuation.  A  cursory  examination  of  waveforms  from  NTS  and  E.  Kazakh 
explosions  shows  that  the  latter  are  much  richer  in  high  frequency  energy.  Some 
attenuation  bias  is  almost  certainly  present,  but  precisely  how  large  is  it? 

APPROACH 

Our  approach  involves  the  following  steps: 

A.  We  analyze  recordings  of  U.S.  and  Soviet  explosions  using  data  from  the 
arrays  EKA,  YKA,  WRA,  GBA  and  NORSAR.  Our  technique  is  to  average  over 
array  elements  and  over  many  similar  events  at  a  test  site.  The  result 
is  smooth  and  consistent  spectra  that  display  the  combined  effect  of 
the  source  and  attenuation  (t*(f)). 

B.  Using  these  spectra  we  develop  t  (f)  models  that  explain  the  spectral 
character  from  below  0.5  Hz  to  as  high  as  8  Hz  (for  E.  Kazakh). 

C.  We  explore  tradeoffs  between  source  and  attenuation  models. 

D.  The  m^  bias  is  estimated  with  synthetic  seismograms  that  provide  a  good 
match  to  the  data  in  both  the  time  and  frequency  domain. 

E.  Our  results  are  compared  to  results  from  other  studies  and  the  conclu¬ 
sions  summarized  within  a  broad  context. 

This  approach  differs  from  that  of  previous  studies  by  using  spectra  that  are 
more  consistent  and  simpler  to  interpret.  Also,  our  interpretation  is  based  on 
synthetic  seismograms  that  force  a  specific  statement  of  all  assumptions.  Thus, 
there  is  no  way  to  avoid  confrontation  with  the  uncertainties.  Finally,  we  have 
tried  to  be  reasonably  complete,  synthesizing  much  of  the  relevant  evidence  in 
one  study. 


RESULTS 


The  array  spectrum  for  each  event  is  obtained  by  averaging  spectra  from  the 
individual  array  elements.  For  each  element  we  compute  the  energy  density 
spectra  for  short  time  windows  isolating  the  P  wave,  then  subtract  the  power 
spectrum  for  a  noise  sample  taken  before  P.  These  "noise-corrected"  energy 
density  spectra  are  then  averaged  over  all  elements  of  the  array  and  converted  to 
Fourier  spectra. 

The  array  spectra  variations  are  consistent  with  yield- dependent  corner  frequency 
variations.  Correcting  for  the  source,  we  can  deduce  the  attenuation.  We  assume 
that  the  NTS  explosions  can  be  represented  by  the  Mueller  and  Murphy  (1971) 
tuff/rhyolite  source  model  while  E.  Kazakh  explosions  are  represented  by  the 
Mueller/Murphy  granite  model.  We  make  "source  corrections"  to  each  event 
spectrum  with  these  assumptions  and  average  the  "source-corrected"  spectra  for 
groups  of  similar  events.  The  results  are  shown  in  Figure  2  for  E.  Kazakh  and  in 
Figure  3  for  NTS.  The  m^  bias  depends  on  attenuation  below  about  2  Hz,  but  this 
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Figure  2  Source  corrected  path  spectra  are  plotted  for  each  data  set  considered. 

In  (a)  we  compare  the  path  spectra  from  different  portions  of  the  test 
site  at  each  UKAEA  station  (at  YKA  there  are  too  few  data  to  bisect  the 
Shagan  River  site) .  In  (b)  the  mean  spectrum  for  the  Degelen  paths  is 
shown  with  dotted  lines  indicating  the  sample  standard  deviation.  At 
the  bottom  of  (b)  la  shown  a  comparison  of  the  mean  spectra  for  the  paths 
from  that  particular  source  area  to  each  station.  (Reproduced  from  Bache 
et  al.,  1986). 


Figure  3  The  path-average  spectra  for  Piledriver  and  five  event  populations  are 
compared.  The  number  in  parentheses  indicates  the  number  of  events  in 
each  population.  In  each  of  these  comparisons  the  spectra  have  been 
aligned  so  they  coincide  at  high  frequencies. 
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depends  on  our  assumptions  about  the  source.  This  is  indicated  in  Figure  4, 
where  the  source-corrected  path  spectra  are  shown  for  several  choices  for  the 
source  model. 

For  any  reasonable  source  assumptions,  it  is  clear  that  frequency -dependent 

attenuation  is  required  to  match  the  data.  We  use  a  two-parameter  absorption 

band  model.  The  parameters  are  t*  (the  long  period  t*)  and  rm,  which  controls 
the  frequency- dependence .  A  variety  of  models  that  fit  the  data  are  shown  in 
Figure  5.  All  give  time  domain  waveforms  that  match  the  data  equally  well.  The 
for  various  synthetic  seismograms  are  plotted  in  Figure  6  with  the  "m^  bias" 
due  to  attenuation  effects  only.  There  can  be  additional  bias  due  to  source - 
coupling  differences.  For  example,  the  i‘;fferences  between  Hue  Her -Murphy  tuff 
and  granite  sources  introduces  an  add?  onal  bias  of  0.2  m^.  However,  there  is 
little  or  no  evidence  that  source -coupling  differences  this  large  occur  at  NTS, 
much  less  at  E.  Kazakh. 

The  "medium  bias"  results  are  most  consistent  with  conventional  ideas  for  source 
spectra  and  the  attenuation  differences  at  long  period  between  NTS  and  shield¬ 
like  regions.  The  "small  bias"  can  occur  if  the  E.  Kazakh  explosions  have 

substantially  more  high  frequency  energy  than  the  Mueller -Murphy  granite  source, 
or  if  the  apparent  attenuation  differences  near  1  Hz  reflect  rm  differences  with 
At*  rather  small. 


Figure  4  The  effect  of  several  al¬ 
ternative  source  correc¬ 
tions  are  shown  for  the 
GBA  S.W.  Shagan  spectra. 
The  amplitudes  are 
normalized  to  be  identical 
at  high  frequency.  The 
M-Ml  is  the  standard 
Mueller-Murphy  granite 
source.  The  VS-B  is  the 
von  Seggern  and  Blandford 
(1972)  granite  source  with 
the  same  assumptions  (m^ 
6.0  corresponds  to  150Kt 
and  mb*  0.9  lor  W).  The 
M-M2  and  M-M3  sources  are 
nearly  the  same  as  M-Ml, 
except  that  the  frequen¬ 
cies  have  been  shifted 
down  (M-M2 )  or  up  (M-M3 ) 
by  a  factor  of  1.6. 
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Figure  5  Source-corrected  EKA  Pahute  Mesa  and  GBA  Shagan  Ri^er  spectra  are 

compared  with  the  spectral  effect  of  the  indicated  attenuation  model. 
For  the  GBA  Shagan  data  attenuation  models  are  shown  for  two  choices 
for  the  source  model.  Also  shewn  in  three  of  the  plots  is  the  spectrum 
of  a  synthetic  seismogram  (processed  the  same  way  as  the  data)  with  the 
indicated  P-pP  lag  time. 
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Estimates  of  Bias  Due  to  Attenuation  Differences 
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Figure  6  Synthetic  are  plotted  versus  yield.  The  line  is  a  fit  to  n^-yield 

data  for  NTS  explosions  below  the  water  table  (Murphy,  1981).  The 
differences  are  tabulated  at  the  top. 
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There  are  many  other  lines  of  evidence  that  can  be  pursued  to  help  decide  which 
bias  estimate  should  be  preferred.  Key  elements  include  the  following: 

o  Our  results  are  entirely  consistent  with  those  obtained  by  Murphy 
(1985)  with  network- average  spectra.  However,  our  spectra  are  broader- 
band,  forcing  us  to  include  frequency- dependence  in  the  attenuation, 
and  we  consider  more  variations  in  the  possible  source  spectrum. 

o  There  are  large  uncertainties  in  the  source  spectrum  for  explosions  in 
hard  rock.  The  data  from  U.S.  and  French  Sahara  explosions  are  not 
adequate  to  exclude  the  Mueller -Murphy  granite  source.  However,  they 
are  easier  to  explain  if  the  French  Sahara  source  is  smaller  than  the 
U.S.  granite  source,  and  if  the  U.S.  granite  source  has  larger  over¬ 
shoot  (smaller  4'<l))  than  conventional  models. 

o  The  best  available  analysis  of  the  surface  wave  data  (Stevens,  1986) 
concludes  that  a  source  spectrum  with  large  overshoot  is  more  consis¬ 
tent  with  the  E.  Kazakh  data. 

o  Direct  measurement  of  m^  bias  between  NTS  and  sites  on  the  Canadian 
Shield  should  give  an  upper  bound  for  the  NTS  -  E.  Kazakh  bias.  The 
SDCS  experiment  (Der  et  al.,  1981)  provides  a  reciprocal  estimate  for 
this.  These  results  can  only  be  used  to  support  m^  bias  greater  than 
0.2  by  constructing  ad  hoc  arguments  for  discarding  the  measured  m^ 
bias  in  favor  of  "inferred"  or  "corrected"  values. 

o  Taken  at  face  value,  the  Shagan  River  cratering  explosion  indicates  an 
that  150Kt  corresponds  to  m^,  5.92  (total  bias  of  about  0.25).  There  is 
no  compelling  evidence  to  suggest  that  this  should  be  corrected  upward 
by  more  than  0.1  for  the  difference  between  cratering  and  contained 
explosions  (Day  et  al.,  1986).  We  cannot  use  any  one  event  to  estimate 
bias  with  much  confidence,  but  note  that  this  event  cannot  be  used  as 
evidence  to  support  large  m^  bias  estimates. 

We  also  computed  spectra  for  explosions  at  the  French  Sahara  and  Mururoa  test 
sites.  The  source  spectrum  assumptions  again  control  conclusions  about  at¬ 
tenuation.  The  simplest  interpretation  is  that  attenuation  at  the  Sahara  site  is 
intermediate  between  that  for  NTS  and  Shagan  River.  The  attenuation  for  Mururoa 
appears  to  be  greater  (higher  t*)  than  for  NTS. 


CONCLUSIONS 

We  conclude  that  there  is  no  way  of  conclusively  proving  either  hypothesis  listed 
in  the  Introduction  without  more  information  about  the  source  coupling  of  the 
Shagan  River  explosions.  The  second  hypothesis  (yields  over  150Kt  at  Shagan 
River)  easily  explains  the  data  if  the  Shagan  River  source  spectra  are  different 
(more  overshoot,  higher  corner  frequency)  from  conventional  source  models  for 
explosions  in  hard  rock.  But  these  conventional  models  have  not  been  shown  to 
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match  the  hard  rock  data,  and  these  data  may  be  easier  to  explain  with  lower 
amplitude,  higher  overshoot  source  spectra  of  the  kind  that  lead  to  the  lower  mb 
bias  estimates. 

It  is  possible  to  start  with  either  one  of  the  two  hypotheses  (or  many  others) 
and  work  backward  to  fit  all  the  data  within  a  framework  consistent  with  that 
hypothesis.  When  data  do  not  fit,  plausible  ad  hoc  corrections  can  be  introduced 
to  explain  the  discrepancy.  The  data  are  often  ambiguous,  and  there  are  many 
things  (e.g.,  source  coupling  and  frequency-dependent  attenuation)  we  do  not 
understand  very  well.  Taking  all  factors  into  consideration,  the  low  mb  bias 
hypothesis  seems  more  consistent  with  all  the  data  in  the  sense  that  fewer  ad  hoc 
arguments  are  required  to  explain  away  data  inconsistent  with  this  hypothesis. 
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ABSTRACT 

It  has  long  been  realized  that  oceanic  Pn  and  Sn  phases  retain  extremely  high  frequen¬ 
cies  even  after  propagation  to  teleseismic  distances.  Estimates  of  the  effective  quality 
factor,  Q  ,  indicate  a  nearly  elastic  rheology  for  the  oceanic  lithosphere.  Data  collected 
during  the  1983  Ngendei  Seismic  Experiment  in  the  southwest  Pacific  are  consistent 
with  these  observations  at  high  frequencies,  but  deviate  substantially  at  low  frequencies. 
Spectra  of  Pn  data,  when  plotted  on  a  logarithmic  scale,  display  an  obvious  break  in 
slope  near  a  frequency  of  1-2  Hz.  When  Q  is  estimated  from  the  spectral  slope  of  these 
events,  a  much  lower  value  is  obtained  at  frequencies  less  than  1-2  Hz  than  at  frequen¬ 
cies  in  excess  of  5-10  Hz.  While  an  array  of  instruments  was  not  available  to  eliminate 
the  source  spectra  by  the  spectral  ratio  method ,  any  of  the  standard  source  spectra 
models  would  yield  this  result.  The  choice  of  the  high  frequency  fall-off  rate  and  the 
corner  frequency  strongly  influence  the  values  obtained  for  Q ,  but  not  the  observation 
that  Q  increases  with  frequency.  Similar  observations  for  travel  paths  in  the  northwest 
Pacific  recorded  by  the  Wake  Island  hydrophone  array  have  also  been  reported.  In 
order  to  synthesize  more  realistic  oceanic  Pn  and  Sn  wave  trains,  it  now  appears  that 
Q  should  increase  with  frequency.  Wave  number  integration  is  conducted  in  the  fre¬ 
quency  domain  and  attenuation  is  introduced  by  allowing  material  velocities  to  be  com¬ 
plex.  We  have  adopted  an  absorption  band  rheology  with  a  distribution  of  relaxation 
times  to  implement  a  frequency  dependence  of  Q  into  the  wave  number  integration 
algorithm.  An  absorption  band  consistent  with  published  Q  values  and  the  Ngendei 
spectra  was  used.  Synthetic  Pn  and  Sn  phases  were  generated  to  a  frequency  of  6.4  Hz 
on  a  Cray  X-MP  computer.  The  results  of  the  modeling  indicate  a  much  better  fit  of 
synthetics  to  data  is  possible  by  allowing  Q  to  be  a  function  of  frequency. 

INTRODUCTION 

Oceanic  Pn  and  Sn  phases  are  commonly  observed  on  ocean  bottom  seismometer 
(OB^)  records  in  the  distance  range  of  5  '  to  30  ’ .  These  are  guided  phases  that 
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propagate  entirely  within  the  oceanic  lithosphere  and  are  characterized  by  an  onset 
velocity  of  8.0  km/s,  a  long  reverberatory  wave  train,  and  the  presence  of  extremely 
high  frequencies  even  for  regional  to  teleseismic  path  lengths.  Frequencies  as  high  as  15 
Hz  for  Pn  and  20  Hz  for  Sn  have  been  reported  at  a  distance  of  3300  km  for  travel 
paths  under  the  western  Pacific  (  Walker ,  et  al .,  1983).  Walker,  et  al.  (1978)  estimated 
the  effective  quality  factors  for  Pn  and  S„  with  travel  paths  in  the  northwest  Pacific  for 
two  separate  earthquakes  yielding  P%  values  of  3700±200  and  8400±1300  and  SH  values 
of  8500±600  and  19,100±3700.  More  recently,  Butler  et  al.  (1985)  report  a  strong  fre¬ 
quency  dependence  of  QPn  and  Q$n  between  2  and  20  Hz  for  travel  paths  in  the 
northwest  Pacific.  They  were  able  to  use  the  spectral  ratio  method  between  pairs  of 
receivers  to  estimate  the  quality  factors  for  Pn  and  Sn .  They  reported  Q  values  of  300 
at  2.5  Hz  and  1500  at  17.5  Hz  for  Pn  and  for  S„  ,  500  at  2.5  Hz  and  3000  at  22  Hz. 

SYNTHETIC  Pn  AND  Sn  PHASES 

Wave  number  integration  was  used  to  generate  complete  synthetic  seismograms  for 
an  oceanic  lithosphere  model  with  a  constant  Q  rheology.  The  velocity  structure  was 
chosen  on  the  basis  of  a  refraction  study  conducted  during  the  1983  Ngendei  Seismic 
Experiment  in  the  southwest  Pacific  (23  ‘49 'S,  165*32'  W).  The  Ngendei  site  was  the 
location  of  hole  595B  of  Deep  Sea  Drilling  Project  (DSDP),  leg  91.  The  depth  of  the  oce¬ 
anic  water  column  at  the  site  is  5.55  km  and,  from  core  sample  measurements,  the  sedi¬ 
ment  thickness  is  70  m  with  a  compressional  velocity  of  1.6  km/s.  The  model  is  listed  in 
Table  1  and  a  single  source  depth  of  15  km  was  used. 

Figure  1  illustrates  vertical  and  radial  displacement  at  the  sea  floor  at  a  range  of 
1000  km  for  a  sub-Moho  thrust  fault.  The  fault  used  had  a  dip  of  25  ' ,  a  rake  of  270  * , 
and  an  azimuth  relative  to  the  receiver  of  60  ’ .  The  synthetic  radial  component  oceanic 
Pn  and  Sn  phases  very  closely  resemble  data  collected  at  the  Ngendei  site.  The  substan¬ 
tial  coda  of  these  phases  is  modeled  as  a  sum  of  leaky  organ-pipe  modes  in  the  oceanic 
water  column  and  sediment  layer.  Shear  wave  resonance  is  particularly  prominent  on 
the  horizontal  component  due  to  the  steep  incidence  of  the  rays  at  the  receiver.  The 
’pulse-like’  character  of  the  water  reverberations  in  the  synthetic  vertical  component  Pn 
wave  train  does  not  ,  however,  emulate  the  data.  The  vertical  component  coda  results 
predominantly  from  compressional  wave  reverberation  in  the  oceanic  water  column  and 
sediment  layer.  Due  to  the  low  impedance  contrast  at  the  ocean-sediment  interface  for 
compressional  waves,  the  sediment  modes  are  extremely  leaky  and  die  out  long  before 
the  arrival  of  the  next  water  multiple. 

Figure  2a  is  the  synthetic  vertical  displacement  of  Figure  la  after  including  the 
OBS  instrument  response  and  a  source  spectrum  with  /  ~2  amplitude  decay  and  corner 


203 


frequency  of  0.6  Hz.  Figure  2b  is  the  vertical  component  Pn  spectrum.  The  sharp 
equally  spaced  peaks  are  due  to  resonance  within  the  oceanic  water  column.  Figure  2c  is 
an  example  of  vertical  component  data  collected  at  the  Ngendei  site  and  Figure  2d  is  the 
first  12.8  Hz  of  its  power  spectrum.  The  obvious  difference  between  the  two  spectra  is 
the  relative  enrichment  of  the  power  at  frequencies  between  0-2  Hz  of  the  data.  The  lack 
of  this  low  frequency  energy  gives  rise  to  the  ’pulse-like’  appearance  of  the  synthetic  Pn 
wave  train. 


ABSORPTION  BAND  Q 

One  way  to  enhance  the  power  level  of  the  low  frequencies  relative  to  higher  fre¬ 
quencies  is  to  lower  the  Q  .  This,  however,  inhibits  the  propagation  of  extremely  high 
frequencies  to  teleseismic  distances  required  by  the  data.  As  an  alternative,  Q  could 
increase  with  frequency  thereby  enhancing  the  power  at  frequencies  between  0-2  Hz 
relative  to  frequencies  between  3-5  Hz  without  severely  attenuating  frequencies  in  excess 
of  15  Hz.  To  approximate  the  true  frequency  dependence  of  Q  in  the  oceanic  litho¬ 
sphere,  we  have  adopted  an  absorption  band  rheology  with  a  distribution  of  relaxation 
times.  This  is  incorporated  into  the  wave  number  integration  algorithm  ( Lundquiat  and 
Cormier,  1980)  by  allowing  the  complex  material  velocity  to  be  of  the  form 
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where  V,  is  the  ’elastic’  velocity,  Qm  is  the  minimum  Q,  and  rx  and  r2  are  the  low 
and  high  frequency  relaxation  times,  respectively.  Care  must  be  taken  in  applying  this 
equation  in  that  the  natural  log  term  perturbs  the  real  part  of  the  complex  velocity  and 
can  cause  an  unrealistic  deviation  in  travel  times.  In  the  limit  as  this  term  approaches 
zero,  the  familiar  expression  for  Q  as  a  function  of  frequency  for  an  absorption  band 
rheology  is  obtained 
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As  an  example  of  the  effect  of  an  absorption  band  Q  on  the  spectra  we  used 
values  for  compressional  waves  of  Qm  =200  and  r2=0.05  s  resulting  in  a  linear  frequency 
dependence  of  Q  for  frequencies  greater  than  about  3.0  Hz.  The  low  frequency  relaxa¬ 
tion  time,  ru  was  arbitrarily  set  to  100s,  well  out  of  the  frequency  range  of  interest.  Fig¬ 
ure  3a  displays  curves  of  the  instrument  response  times  a  source  spectrum  with  f  ~2 

-rtf 

amplitude  decay  and  corner  frequency  of  0.6  Hz  multiplied  by  e  for  a  range  of 
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1000  km,  a  group  velocity  of  8  km/s,  and  varying  Q  values.  The  dashed  curves 
represent  constant  Q  values  and  for  the  solid  curve  the  absorption  band  Q  was  used. 
Figure  36  displays  the  same  curves  to  a  frequency  of  64  H2.  This  illustrates  that  to 
enhance  the  low  frequency  energy  of  the  synthetics  without  depleteing  the  very  high  fre¬ 
quencies,  Q  must  increase  with  frequency. 

Synthetic  vertical  and  radial  displacement  for  the  Ngendei  model  with  an  absorp¬ 
tion  band  rheology  in  the  sub-Moho  portion  of  the  lithosphere  are  shown  in  Figure  4. 
While  the  modeling  is  not  yet  complete,  the  agreement  of  the  synthetic  vertical  com¬ 
ponent  PH  wave  train  with  the  data  is  greatly  enchanced  by  allowing  Q  to  increase 
with  frequency. 
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Table  1.  Ngendei  Model 


Layer 

Vf  ,  km/s 

V, ,  km/s 

p,  Mg/m3 

Qa 

h,  km 

1 

1.500 

0.000 

1.0000 

50000 

— 

5.55 

2 

1.600 

0.116 

1.0096 

200 

100 

0.07 

3 
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0.25 
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5 
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6 
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0.20 

7 

6.200 

3.419 

2.9400 
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0.40 

8 

6.500 

3.585 

3.0000 
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0.70 

9 

6.800 

3.750 

3.0600 

900 

450 

4.63 

10 

8.200 
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3.3400 

2000 
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11 
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3.3450 

2000 
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12 
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3.3500 

2000 
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3.3550 

2000 
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Figure  1.  Synthetic  Pn/  Sn  phases  computed  for  the  Ngendei  structure,  (a)  Vertical 
displacement  and  (6)  radial  displacement  for  a  thrust  fault  source  at  15.0  km  depth  at 
an  epicentral  range  of  1000  km.  Velocities  of  8.1  and  4.65  km/s  are  indicated. 
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Figure  2.  (a)  Synthetic  vertical  component  Pn  phase  of  Figure  la  after  including  the 
OBS  instrument  response  and  a  source  spectrum  with  /  2  amplitude  decay  and  corner 
frequency  of  0.6  Hz.  ( b )  Power  spectrum  of  the  Pn  phase  in  Figure  2a.  (c)  Vertical  com¬ 
ponent  Pn  data  taken  at  the  Ngendei  site.  (</)  First  12.8  Hz  of  the  data  spectrum. 
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Figure  3.  Curves  of  instrument  response  times  a  source  spectrum  with  /  ~2  amplitude 

ZllL 

decay  and  corner  frequency  of  0.6  Hz  multiplied  by  e  for  a  range  of  1000  km,  a 
group  velocity  of  8  km/s,  and  varying  Q  values.  The  dashed  curves  are  for  constant  Q 
and  the  solid  curve  represents  the  result  for  an  absorption  band  Q  model,  (a)  First  12.8 
Hz  and  (4)  64.0  Hz  of  the  spectrum. 
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Figure  4.  Synthetic  Pn/  Sn  phases  computed  for  the  Ngendei  structure  with  an  absorp¬ 
tion  band  Q  adopted  in  the  mantle  portion  of  the  lithosphere,  (a)  Vertical  displace¬ 
ment  and  ( b )  radial  displacement  for  a  thrust  fault  source  at  15.0  km  depth  at  an  epi- 
central  range  of  1000  km.  Velocities  of  8.1  and  4.65  km/s  are  indicated. 
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SURFACE  WAVE  SYNTHESIS  FOR  LATERALLY  HETEROGENOUS  SPHERICAL  EARTH  MODELS 


Keiiti  Aki,  Ta-liang  Teng,  John  Faulkner 
Department  of  Geological  Sciences,  University  of  Southern  California 

INTRODUCTION 

The  subject  matter  of  our  research  deals  with  the  improvement  on  the  absolute 
accuracy  concerning  the  Ms:yield  determination.  At  large  epicentral  distances,  surface 
waves  are  the  most  prominent  phase  in  a  seismogram.  If  they  are  to  be  effectively 
used  for  treaty  monitoring  purposes,  uniform  determination  of  Ms  to  an  accuracy  of 
0.1  magnitude  unit  becomes  important.  This  translates  to  an  absolute  determination 
of  surface-wave  amplitudes  (in  time  or  frequency  domain)  to  an  accuracy  of  better 
than  25%.  Factors  entering  into  any  amplitude  measurement  must  include  the  source 
effect,  the  propagation  effect,  and  the  recording  site  effect;  each  of  these  effects 
makes  a  significant  contribution  to  data  scatter.  All  of  these  effects  are  currently 
under  intensive  study  and  these  research  results  should  lead  to  a  more  accurate  yield 
determination  as  well  as  discrimination. 

Over  the  past  12-month  contract  period,  we  have  accomplished  the  following  tasks 
with  results  that  have  considerably  improved  our  understanding  on  the  propagation 
of  surface  waves  in  a  laterally  heterogeneous  earth.  The  nature  of  focusing  and 
defocusing  of  the  surface  wavefront  has  important  bearings  on  the  basic  analysis 
of  dispersion  studies  and  therefore  the  structural  inversion.  Its  implication  on 
surface-wave  amplitude  variations,  of  course,  is  the  main  thrust  of  this  contract 
research . 

GAUSS IAN-BEAM  SYNTHESIS  OF  RAYLEIGH  WAVES 

For  the  wave-form  synthesis  of  Rayleigh  waves  in  a  laterally  heterogeneous  Earth, 
we  use  the  Gaussian  beam  method  originally  developed  for  body  waves  by  Cerveny  et 
al .  (  1982)  and  modified  for  surface  waves  by  Yomogida  and  Aki  (  1985). 

The  following  is  a  brief  summary  of  how  to  make  a  synthetic  seismogram  of  surface 
waves  by  the  Gaussian  beam  approach.  First,  we  specify  the  phase  velocity  at  mesh 
points  and  apply  the  transformation  corresponding  to  the  Mercator  projection  (Jobert 
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and  Jobert,  1983).  Second,  we  apply  the  bi-cubic  spline  interpolation  of  phase  velocity 
so  that  we  can  calculate  its  first  and  second  derivative  at  any  point.  Third,  we 
solve  the  r  y  tracing  equation  and  shoot  ray  paths  from  a  given  source  to  all  azimuths 
at  2°  or  finer  intervals.  Fourth,  we  solve  the  dynamic  ray-tracing  equation  and 
determine  the  spreading  and  wave-front  curvature  for  each  ray.  Fifth,  we  construct 
the  Gaussian  beam  solution  for  each  ray  corresponding  to  a  source  function  of  the 
Gabor  type  and  sum  up  the  contribution  of  each  Gaussian  beam  to  a  particular  station. 
Sixth,  we  determine  the  amplitude  and  arrival  time  of  wave  packets  using  the  group 
velocity  data  and  the  kinetic  energy  integral  calculated  for  the  appropriate  model 
structure . 

The  computer  program  for  the  synthesis  has  been  written  by  K.  Yomogida  at  MIT 
and  was  adapted  to  the  USC  computer  by  J.  Faulkner.  A  number  of  test  cases  have 
been  run  to  make  sure  that  the  code  gives  correct  results.  A  world  map  of  Mercator 
projection  has  also  been  digitized  that  forms  the  base  of  reference  for  the  velocity 
map  and  for  the  nuclear  test  sites  and  recording  stations. 

In  order  to  carry  out  the  above  steps  for  ray  tracing  and  synthetic  seismogram 
computations,  we  first  need  a  world  map  of  phase  velocity  for  each  period  of  interest. 

WORLD  MAP  OF  RAYLEIGH-WAVE  PHASE  VELOCITY 

To  perform  the  ray-tracing  and  Gaussian-beam  synthesis,  phase  velocity,  group 
velocity,  and  surface-wave  energy  integral  data  are  needed  for  the  entire  region 
in  which  waves  propagate.  Recently,  Rosa  (personal  communication)  collected  all 
the  existing  phase  velocity  data  (published  and  unpublished)  for  the  period  range 
from  20  to  100  sec.  For  example,  he  collected  phase  velocity  at  period  of  40  sec. 
for  790  paths. 

Rosa  compared  these  ray  paths  with  the  regionalized  tectonic  map  of  Jordan  (1981), 
and  selected  those  paths  which  have  at  least  70%  of  the  total  path  length  lying  within 
one  of  Jordan's  six  regions,  namely  (A)  young  ocean  (0-25  My),  (B)  intermediate  ocean 
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(25-100  My),  (C)  old  ocean  (  -’100  My),  (P)  Phanerozoic  platform,  (Q)  Phanerozoic 
orogenic  zones,  and  (S)  Precambrian  shields  and  platforms.  He  then  computed  the 
mean  and  the  standard  error  for  each  region  for  periods  20,  30,  40,  50,  60,  70,  80, 

90,  and  98  sec.  We  assigned  these  phase  velocity  values  at  the  center  of  the  5° 
x  5°  area  which  is  classified  to  one  of  the  above  six  regions  by  Jordan  (Fig.  3). 

For  this  preliminary  study,  we  assumed  that  the  group  velocity  is  constant  (4.0 
km)  for  all  regions  and  all  periods.  This  affects  the  arrival  time  of  the  wave  packet 
significantly  but  only  slightly  the  amplitude  of  individual  Gaussian  beam  solutions. 

EXAMPLES  OF  SURFACE-WAVE  RAY  MAPS  AND  THE  CORRESPONDING  SYNTHETIC  RECORDS 

Based  on  Rosa's  global  phase  velocity  data,  we  have  computed  test  cases  for 
a  number  of  surface-wave  ray  maps  and  the  corresponding  Gaussian-beam  synthetic  seismograms. 
We  have  generated  ray  maps  for  periods  from  20  sec  to  90  sec  at  10  sec  increments. 

Figures  1-abc  give  ray  maps  for  30  sec  surface  waves  originated  from  Novaya  Zemlya, 

East  Kazakh,  and  China,  respectively.  It  is  interesting  to  see  that  for  an  assumed 
velocity  model  of  mild  lateral  heterogeneity,  intermediate  period  surface  waves 
can  significantly  alter  their  wave  paths.  Focusing  and  defocusing  regions  as  well 
as  "transmission  corridors"  are  well  displayed  that  shift  themselves  with  the  change 
of  source  locations.  Unfortunately,  because  of  the  nature  of  the  Mercator  projection, 
the  path  lying  near  the  North  Pole  cannot  be  included  here.  Furthermore,  ray  paths 
terminated  at  the  western  boundary  of  the  map  are  not  included  in  the  Gaussian-beam 
synthesis.  Therefore,  the  synthetic  seismograms  as  now  presented  in  Figure  2  have 
errors  in  the  R1  for  eastern  U.S.  stations  (where  Rl  is  missing)  and  in  R2  for  western 
U.S.  stations  (where  R2  is  missing).  This  type  of  error  is  being  corrected  now, 
which  involves  a  coordinate  rotation  of  the  North  Pole  out  of  the  way  to  perhaps 
the  equator  with  an  attendant  mapping  of  the  global  phase  velocity  distribution  onto 
the  new  rotated  coordinates. 

We  are  currently  checking  into  actual  records  to  see  if  the  ray  maps  correctly 
predict  the  observations. 
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Figure  2 

E  Kazakh  (0  =50°  X  =78°) 

(T  =  50  Seconds,  L,  =  1.0,  SQ  =  1.0,  <5.  =2°,Q=1000) 


(spuooasjauiij, 


o 

o 

o 

o 

o 

o 

o  o  o 

o 

O 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

o  o  o 

o 

o 

o 

o 

o 

o 

o 

o 

in 

o 

in 

o 

IT) 

o  m  o 

in 

o 

in 

o 

m 

o 

in 

O) 

CO 

oo 

Pv 

U3 

to  m  m  , 

•4- 

m 

fO 

CN 

CM 

(  Peak  Amplitude  =  0  81E-02  ) 

Station  Number 


AFGL  /  DARPA  Seismic  Research  Symposium 
U.  S.  Air  Force  Academy 
Colorado  Springs  -  6  -  8  May  1986 


Paper  tittle  :  Regional  Distribution  of  apparent  attenuation 

in  the  Center  of  France 


Paper  authors  :  M.  Campillo  -  Radiomana 


Grant  n°  85  0033 


The  set  of  seismic  data  is  the  one  already  used  in  our  previous  studies 
on  Lg  '*'aves  (  18  quakes  with  ML  between  3*2  and  4-8,  within  the  crust 
22  SP  stations  of  the  LDG  network  ). 

The  method  used  here  is  the  Algebric  Reconstruction  Technic  (ART) 
to  build  up  an  image  of  the  Q  factor  spatial  distribution. 

The  set  of  propagation  paths  is  presented  in  Figure  1. 

Each  point  of  the  image  of  Q.  takes  into  account  all  the  paths  within 
a  distance  of  50  km. 

In  order  to  apply  some  confidence  criterium  on  the  quality  of  this  Q 
evaluation  we  define  a  confidence  index  : 

i  c  =  n  source  x  1  rays 

where  : 

n  source  is  the  numbers  of  different  sources  which  generate 

the  rays 

1  rays  the  length  of  the  ray  inside  the  influence  circle  divided  by 
the  influence  radius. 

A  map  of  this  index  is  given  in  Figure  2  ,  where  the  represented  points 
fit  the  test  of  ray  length. 
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As  a  result,  all  the  blue  zones  correspond  to  a  high  confidence  index, 
and  red  zones  might  include  systematic  bias. 

Figure  3  represents  a  rather  simple  regional  distribution.  If  one  just 
keeps  within  high  confidence  limits,  the  main  regions  are  defined  as  : 

-  Limagne  Basin  ,  Forez,  Saone  valley  =  weak  figure  for  Q 

(  between  Q °  =  180  and  Q 0  =  240  ) 

-  Limousin  =  (  between  Q»  =  200  and  Q«  =  260  ) 

-  South  Parisian  basin,  Berry,  Poitou,  North  of  Massif  Central 
(  between  Q°  =  300  and  Qo  =  440  )• 

The  mean  Q°  value  (  always  computed  at  f  =  1  Hz  )  extracted 

from  these  results  is  :  Q>  =  290 

Similar  to  the  one  obtained  by  direct  Lg  measurement  (  Campillo 
and  al  1984  ). 

We  need  to  extend  this  Q°  cartography  to  other  western  Europe  regions 
by  using  other  stations  and  earthquakes. 


Conclusion 


These  results  pointed  out  the  importance  of  graben  zones  for  Lg  waves 
propagation.  Nevertheless,  at  this  stage  of  the  study  we  do  not  possess 
yet  any  separation  tool  between  two  types  of  phenomena  which  could 
create  attenuation 

-  small  scale  diffraction  (  or  scattering  ) 

-  large  scale  geometrical  effect 

A  study  of  frequency  dependancy  for  Q  should  help  us  in  the 
interpretation  of  our  data. 

On  the  other  hand,  developpements  of  numerical  simulation  methods 
encourage  us  to  investigate  more  complex  zones  in  order  to  interprete 
the  quasi  extinction  phenomena  of  Lg  waves  we  have  observed  along 
peculiar  paths. 
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Figure  l 
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